Imported Upstream version 0.1.6~dfsg

diff --git a/ChangeLog b/ChangeLog
index 3bf82a59c9f2f5d0249c9d457022a93f3758d751..c0afc458da745949eaecc674e2adc4c264f54e6a 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,561 @@
+------------------------------------------------------------------------
+r451 | lh3lh3 | 2009-09-02 10:44:48 +0100 (Wed, 02 Sep 2009) | 4 lines
+Changed paths:
+   M /trunk/samtools/bam_md.c
+   M /trunk/samtools/bam_rmdup.c
+   M /trunk/samtools/bam_rmdupse.c
+   M /trunk/samtools/bam_sort.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/samtools.1
+
+ * samtools-0.1.5-34 (r451)
+ * applied the patch by John
+ * improved the help message a little bit
+
+------------------------------------------------------------------------
+r450 | lh3lh3 | 2009-09-02 09:55:55 +0100 (Wed, 02 Sep 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam_color.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-33 (r450)
+ * fixed a bug in bam_color.c (on behalf of Nils Homer)
+
+------------------------------------------------------------------------
+r449 | lh3lh3 | 2009-08-29 20:36:41 +0100 (Sat, 29 Aug 2009) | 4 lines
+Changed paths:
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bam_md.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/misc/samtools.pl
+
+ * samtools-0.1.5-32 (r449)
+ * fillmd: fixed a bug in modifying MD/NM tags
+ * in import, give a warning if the read is aligned but there is no CIGAR.
+
+------------------------------------------------------------------------
+r448 | lh3lh3 | 2009-08-19 09:44:28 +0100 (Wed, 19 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/ChangeLog
+   M /trunk/samtools/bam_pileup.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/misc/wgsim_eval.pl
+
+ * samtools-0.1.5-31 (r448)
+ * fixed an issue when the last CIGAR is I or D
+
+------------------------------------------------------------------------
+r447 | lh3lh3 | 2009-08-17 09:34:57 +0100 (Mon, 17 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam_aux.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-30 (r447)
+ * fixed a bug in bam_aux_get(): 'A' is not checked
+
+------------------------------------------------------------------------
+r446 | lh3lh3 | 2009-08-17 09:33:17 +0100 (Mon, 17 Aug 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/bam_aux.c
+   M /trunk/samtools/bamtk.c
+
+ * 
+
+------------------------------------------------------------------------
+r444 | lh3lh3 | 2009-08-11 10:02:36 +0100 (Tue, 11 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam_sort.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-28 (r444)
+ * bug in "merge -n"
+
+------------------------------------------------------------------------
+r443 | lh3lh3 | 2009-08-11 09:29:11 +0100 (Tue, 11 Aug 2009) | 4 lines
+Changed paths:
+   M /trunk/samtools/bam.c
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-27 (r443)
+ * SEQ and QUAL can be "*"
+ * parse CIGAR "=" and "X" as "M"
+
+------------------------------------------------------------------------
+r442 | lh3lh3 | 2009-08-07 21:56:38 +0100 (Fri, 07 Aug 2009) | 4 lines
+Changed paths:
+   M /trunk/samtools/bam_pileup.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/misc/md5.c
+   M /trunk/samtools/misc/md5.h
+   M /trunk/samtools/misc/md5fa.c
+
+ * samtools-0.1.5-26 (r442)
+ * replace RSA Inc md5.* with ones under permissive lincense
+ * fixed a bug in detecting unsorted bam in pileup
+
+------------------------------------------------------------------------
+r441 | bhandsaker | 2009-08-05 14:41:28 +0100 (Wed, 05 Aug 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/bgzf.c
+   M /trunk/samtools/bgzf.h
+   M /trunk/samtools/bgzip.c
+
+Change copyright notices now that MIT has approved open source distribution.
+
+------------------------------------------------------------------------
+r440 | lh3lh3 | 2009-08-05 10:44:24 +0100 (Wed, 05 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam_stat.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-25 (r436)
+ * in flagstats, do not report singletons if both ends are unmapped
+
+------------------------------------------------------------------------
+r439 | lh3lh3 | 2009-08-04 22:16:51 +0100 (Tue, 04 Aug 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/misc/maq2sam.c
+
+fixed a SERIOUS bug in setting 0x20 flag
+
+------------------------------------------------------------------------
+r438 | lh3lh3 | 2009-08-04 21:50:43 +0100 (Tue, 04 Aug 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/misc/samtools.pl
+
+fixed two minor bugs (suggested by Tim M Storm)
+
+------------------------------------------------------------------------
+r437 | lh3lh3 | 2009-08-04 09:13:24 +0100 (Tue, 04 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/misc/samtools.pl
+   M /trunk/samtools/sam_view.c
+
+ * samtools-0.1.5-24 (r435)
+ * fixed a typo
+
+------------------------------------------------------------------------
+r434 | lh3lh3 | 2009-08-03 10:40:42 +0100 (Mon, 03 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam_tview.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-23 (r434)
+ * in tview, press 'r' to show read names rather than sequences
+
+------------------------------------------------------------------------
+r433 | lh3lh3 | 2009-08-02 19:13:35 +0100 (Sun, 02 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/knetfile.c
+
+ * tried to fixed the buggy FTP random access in Windows. FAILED.
+ * anyway, MinGW seems to have problem with "%lld".
+
+------------------------------------------------------------------------
+r432 | lh3lh3 | 2009-08-02 00:32:07 +0100 (Sun, 02 Aug 2009) | 5 lines
+Changed paths:
+   M /trunk/samtools/Makefile.mingw
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/faidx.c
+   M /trunk/samtools/razf.c
+   A /trunk/samtools/win32/libcurses.a
+   A /trunk/samtools/win32/xcurses.h
+
+ * samtools-0.1.5-22 (r432)
+ * faidx: fixed compitability issue with _WIN32
+ * razf: fixed potential compitability issue with _WIN32
+ * PDCurses support in Windows
+
+------------------------------------------------------------------------
+r431 | lh3lh3 | 2009-08-01 23:34:54 +0100 (Sat, 01 Aug 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/win32/libz.a
+
+replace the GnuWin32 version of libz.a with my own build with MinGW.
+
+------------------------------------------------------------------------
+r430 | lh3lh3 | 2009-08-01 23:21:07 +0100 (Sat, 01 Aug 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/knetfile.c
+
+add comments
+
+------------------------------------------------------------------------
+r429 | lh3lh3 | 2009-08-01 22:41:19 +0100 (Sat, 01 Aug 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/Makefile.mingw
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/knetfile.c
+   M /trunk/samtools/knetfile.h
+
+ * samtools-0.1.5-21 (r428)
+ * knetfile.c is now compatible with mingw-winsock
+
+------------------------------------------------------------------------
+r428 | lh3lh3 | 2009-08-01 00:39:07 +0100 (Sat, 01 Aug 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/Makefile.mingw
+
+simplify MinGW Makefile
+
+------------------------------------------------------------------------
+r427 | lh3lh3 | 2009-08-01 00:30:54 +0100 (Sat, 01 Aug 2009) | 5 lines
+Changed paths:
+   A /trunk/samtools/Makefile.mingw
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bamtk.c
+   A /trunk/samtools/win32
+   A /trunk/samtools/win32/libz.a
+   A /trunk/samtools/win32/zconf.h
+   A /trunk/samtools/win32/zlib.h
+
+ * samtools-0.1.5-20 (r427)
+ * MinGW support. At least SAM<->BAM conversion is working. Other
+   functionality are not tested at the moment.
+ * zlib headers and Windows version of libz.a are included in win32/
+
+------------------------------------------------------------------------
+r426 | lh3lh3 | 2009-07-31 23:32:09 +0100 (Fri, 31 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/sam_view.c
+
+ * samtools-0.1.5-19 (r426)
+ * fixed a bug caused by recent modifications. Sorry.
+
+------------------------------------------------------------------------
+r425 | lh3lh3 | 2009-07-31 23:23:51 +0100 (Fri, 31 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/bgzf.c
+
+compatible with Windows binary files
+
+------------------------------------------------------------------------
+r424 | lh3lh3 | 2009-07-31 10:19:59 +0100 (Fri, 31 Jul 2009) | 5 lines
+Changed paths:
+   M /trunk/samtools/bam_maqcns.c
+   M /trunk/samtools/bam_maqcns.h
+   M /trunk/samtools/bam_plcmd.c
+   M /trunk/samtools/bam_tview.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/misc/samtools.pl
+
+ * samtools-0.1.5-18 (r423)
+ * output additional information in pileup indel lines, for the purepose
+   of debugging at the moment
+ * in tview, optionally allow to treat reference skip as deletion
+
+------------------------------------------------------------------------
+r423 | lh3lh3 | 2009-07-30 22:00:36 +0100 (Thu, 30 Jul 2009) | 2 lines
+Changed paths:
+   A /trunk/samtools/misc/psl2sam.pl
+
+convert BLAT psl to SAM.
+
+------------------------------------------------------------------------
+r422 | lh3lh3 | 2009-07-30 11:24:39 +0100 (Thu, 30 Jul 2009) | 6 lines
+Changed paths:
+   M /trunk/samtools/ChangeLog
+   M /trunk/samtools/bam.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/bgzf.c
+   M /trunk/samtools/bgzf.h
+   M /trunk/samtools/knetfile.c
+   M /trunk/samtools/sam.c
+   M /trunk/samtools/sam_view.c
+
+ * samtools-0.1.5-17 (r422)
+ * fixed a but in knetfile.c when seek type is not SEEK_SET
+ * write an empty BGZF block to every BGZF file
+ * check BGZF EOF marker in bam_header_read()
+ * update ChangeLog
+
+------------------------------------------------------------------------
+r421 | lh3lh3 | 2009-07-30 10:03:39 +0100 (Thu, 30 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bam_plcmd.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/misc/samtools.pl
+   M /trunk/samtools/sam.c
+   M /trunk/samtools/sam.h
+   M /trunk/samtools/sam_view.c
+
+ * samtools-0.1.5-16 (r421)
+ * in view and pileup, load header from FASTA index if the input is SAM.
+
+------------------------------------------------------------------------
+r420 | lh3lh3 | 2009-07-29 09:18:55 +0100 (Wed, 29 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/misc/maq2sam.c
+
+do not set "read 1" if reads are not mapped in the PE mode of maq
+
+------------------------------------------------------------------------
+r419 | lh3lh3 | 2009-07-28 09:52:33 +0100 (Tue, 28 Jul 2009) | 5 lines
+Changed paths:
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/misc/samtools.pl
+   M /trunk/samtools/misc/wgsim_eval.pl
+
+ * samtools-0.1.5-15 (r419)
+ * in sam_open(), return NULL when the file cannot be opened.
+ * make wgsim_eval.pl more robust to imperfect SAM
+ * add "unique" command to samtools.pl
+
+------------------------------------------------------------------------
+r418 | lh3lh3 | 2009-07-24 14:04:19 +0100 (Fri, 24 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/misc/wgsim_eval.pl
+
+skip @header lines in SAM
+
+------------------------------------------------------------------------
+r417 | lh3lh3 | 2009-07-24 12:42:38 +0100 (Fri, 24 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/sam_view.c
+
+ * samtools-0.1.5-14 (r417)
+ * more help in "samtools view" due to the recent changes.
+
+------------------------------------------------------------------------
+r416 | lh3lh3 | 2009-07-24 12:34:30 +0100 (Fri, 24 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam.c
+   M /trunk/samtools/bam.h
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/sam.c
+   M /trunk/samtools/sam.h
+   M /trunk/samtools/sam_view.c
+
+ * samtools-0.1.5-17 (r416)
+ * support import/export SAM with string tags
+
+------------------------------------------------------------------------
+r415 | lh3lh3 | 2009-07-24 11:39:26 +0100 (Fri, 24 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam.c
+   M /trunk/samtools/bam.h
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/sam.c
+   M /trunk/samtools/sam.h
+   M /trunk/samtools/sam_view.c
+
+ * samtools-0.1.5-12 (r415)
+ * FLAG now can be in HEX
+
+------------------------------------------------------------------------
+r414 | lh3lh3 | 2009-07-22 22:03:49 +0100 (Wed, 22 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/kstring.h
+
+fixed a compiling error (thank Ken for fixing it)
+
+------------------------------------------------------------------------
+r412 | lh3lh3 | 2009-07-21 22:19:40 +0100 (Tue, 21 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/kstring.c
+   M /trunk/samtools/kstring.h
+
+Implemented Boyer-Moore search in the kstring library.
+
+------------------------------------------------------------------------
+r409 | lh3lh3 | 2009-07-17 17:10:20 +0100 (Fri, 17 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/bam_index.c
+
+do not include knetfile.h when _USE_KNETFILE is not defined
+
+------------------------------------------------------------------------
+r408 | lh3lh3 | 2009-07-17 15:29:21 +0100 (Fri, 17 Jul 2009) | 5 lines
+Changed paths:
+   M /trunk/samtools/bam_md.c
+   M /trunk/samtools/bam_tview.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/bgzf.c
+
+ * samtools-0.1.5-11 (r408)
+ * force to overwirte existing MD if it is different from the one calculated
+   from fillmd.
+ * bgzf.c: improved the compatibility with Windows headers
+
+------------------------------------------------------------------------
+r407 | lh3lh3 | 2009-07-17 14:46:56 +0100 (Fri, 17 Jul 2009) | 5 lines
+Changed paths:
+   M /trunk/samtools/bam.h
+   M /trunk/samtools/bam_aux.c
+   M /trunk/samtools/bam_md.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/sam.h
+
+ * samtools-0.1.5-10 (r407)
+ * implemented bam_aux_del() to remove a tag
+ * fillmd: generate the NM tag
+ * fillmd: cmd interface improvement
+
+------------------------------------------------------------------------
+r406 | lh3lh3 | 2009-07-16 23:30:40 +0100 (Thu, 16 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/Makefile
+
+Sorry. The old Makefile is for PDCurses...
+
+------------------------------------------------------------------------
+r405 | lh3lh3 | 2009-07-16 23:30:11 +0100 (Thu, 16 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/Makefile
+   M /trunk/samtools/bam_tview.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-9 (r405)
+ * improved the compatibility with PDCurses a little bit
+
+------------------------------------------------------------------------
+r404 | lh3lh3 | 2009-07-16 23:23:52 +0100 (Thu, 16 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/Makefile
+   M /trunk/samtools/bam_tview.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-8 (r404)
+ * compatible with PDCurses
+
+------------------------------------------------------------------------
+r403 | lh3lh3 | 2009-07-16 22:39:39 +0100 (Thu, 16 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/kseq.h
+
+ * samtools-0.1.5-7 (r403)
+ * fixed a bug in kseq.h for binary files (text files are fine)
+
+------------------------------------------------------------------------
+r402 | lh3lh3 | 2009-07-16 11:49:53 +0100 (Thu, 16 Jul 2009) | 4 lines
+Changed paths:
+   M /trunk/samtools/bam_import.c
+   M /trunk/samtools/bam_index.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/bgzf.c
+
+ * samtools-0.1.5-6 (r402)
+ * fixed compiling error when "-D_USE_NETFILE" is not applied
+ * improve portability to MinGW
+
+------------------------------------------------------------------------
+r398 | lh3lh3 | 2009-07-13 10:21:36 +0100 (Mon, 13 Jul 2009) | 3 lines
+Changed paths:
+   A /trunk/bam-lite/bam.h (from /trunk/samtools/bam.h:395)
+   A /trunk/bam-lite/bam_lite.c (from /trunk/samtools/bam_lite.c:395)
+   D /trunk/samtools/bam_lite.c
+
+ * move bam_lite.c to bam-lite
+ * copy bam.h to bam-lite
+
+------------------------------------------------------------------------
+r395 | lh3lh3 | 2009-07-13 10:12:57 +0100 (Mon, 13 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam.h
+   M /trunk/samtools/bam_lite.c
+   M /trunk/samtools/bam_lpileup.c
+   M /trunk/samtools/bam_pileup.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-5 (r395)
+ * added bam_pileup_file() and removed bam_lpileup_file()
+
+------------------------------------------------------------------------
+r394 | lh3lh3 | 2009-07-13 00:35:10 +0100 (Mon, 13 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/knetfile.c
+   M /trunk/samtools/knetfile.h
+
+ * samtools-0.1.5-4 (r394)
+ * http_proxy support in knetfile library (check http_proxy ENV)
+
+------------------------------------------------------------------------
+r393 | lh3lh3 | 2009-07-12 23:57:07 +0100 (Sun, 12 Jul 2009) | 5 lines
+Changed paths:
+   M /trunk/samtools/bam_index.c
+   M /trunk/samtools/bam_tview.c
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/knetfile.c
+   M /trunk/samtools/knetfile.h
+
+ * samtools-0.1.5-3 (r393)
+ * knetfile now supports HTTP (no proxy at the moment)
+ * fixed a potential issue in knetfile on opening ordinary file, although I have
+   not seen the sideeffect so far.
+
+------------------------------------------------------------------------
+r392 | lh3lh3 | 2009-07-12 18:50:55 +0100 (Sun, 12 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/samtools.1
+
+Remove the warning in tview
+
+------------------------------------------------------------------------
+r391 | lh3lh3 | 2009-07-12 18:42:43 +0100 (Sun, 12 Jul 2009) | 3 lines
+Changed paths:
+   M /trunk/samtools/bam_tview.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-2 (r391)
+ * do not show a blank screen when no reads mapped
+
+------------------------------------------------------------------------
+r390 | lh3lh3 | 2009-07-09 14:01:42 +0100 (Thu, 09 Jul 2009) | 4 lines
+Changed paths:
+   M /trunk/samtools/bam.h
+   A /trunk/samtools/bam_lite.c
+   M /trunk/samtools/bamtk.c
+
+ * samtools-0.1.5-1 (r390)
+ * removed useless _IOLIB in bam.h. This should cause no change at all.
+ * added bam_lite.c for light-weight BAM reading
+
+------------------------------------------------------------------------
+r385 | lh3lh3 | 2009-07-07 16:53:29 +0100 (Tue, 07 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/knetfile.c
+
+Release samtools-0.1.5c (fixed a bug in piping)
+
+------------------------------------------------------------------------
+r383 | lh3lh3 | 2009-07-07 11:39:55 +0100 (Tue, 07 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/sam.c
+
+Release samtools-0.1.5b (BUG! so embarrassing!)
+
+------------------------------------------------------------------------
+r381 | lh3lh3 | 2009-07-07 11:20:06 +0100 (Tue, 07 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/Makefile
+   M /trunk/samtools/bam.h
+   M /trunk/samtools/bam_aux.c
+   M /trunk/samtools/bamtk.c
+
+Release samtools-0.1.5a (for compatibility with Bio::DB::Sam)
+
+------------------------------------------------------------------------
+r373 | lh3lh3 | 2009-07-07 10:26:57 +0100 (Tue, 07 Jul 2009) | 2 lines
+Changed paths:
+   M /trunk/samtools/ChangeLog
+   M /trunk/samtools/NEWS
+   M /trunk/samtools/bamtk.c
+   M /trunk/samtools/samtools.1
+
+Release samtools-0.1.5
+

 ------------------------------------------------------------------------
 r372 | lh3lh3 | 2009-07-07 09:49:27 +0100 (Tue, 07 Jul 2009) | 3 lines
 Changed paths:
 ------------------------------------------------------------------------
 r372 | lh3lh3 | 2009-07-07 09:49:27 +0100 (Tue, 07 Jul 2009) | 3 lines
 Changed paths:

diff --git a/Makefile b/Makefile
index 7bb44694e0ab26f8266bdadddb25075159b8b365..450b3abbb83876089fbf1da30df600431a272a6e 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,8 +1,6 @@
 CC=                    gcc
 CC=                    gcc

-CXX=           g++

 CFLAGS=                -g -Wall -O2 #-m64 #-arch ppc
 CFLAGS=                -g -Wall -O2 #-m64 #-arch ppc

-CXXFLAGS=      $(CFLAGS)
-DFLAGS=                -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE #-D_NO_CURSES
+DFLAGS=                -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE -D_CURSES_LIB=1

 LOBJS=         bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o \
                        bam_pileup.o bam_lpileup.o bam_md.o glf.o razf.o faidx.o knetfile.o     \
                        bam_sort.o
 LOBJS=         bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o \
                        bam_pileup.o bam_lpileup.o bam_md.o glf.o razf.o faidx.o knetfile.o     \
                        bam_sort.o


@@ -10,9 +8,10 @@ AOBJS=                bam_tview.o bam_maqcns.o bam_plcmd.o sam_view.o \
                        bam_rmdup.o bam_rmdupse.o bam_mate.o bam_stat.o bam_color.o     \
                        bamtk.o
 PROG=          samtools
                        bam_rmdup.o bam_rmdupse.o bam_mate.o bam_stat.o bam_color.o     \
                        bamtk.o
 PROG=          samtools

-INCLUDES=      
+INCLUDES=

 SUBDIRS=       . misc
 SUBDIRS=       . misc

-LIBPATH=       
+LIBPATH=
+LIBCURSES=     -lcurses # -lXCurses

 
 .SUFFIXES:.c .o
 
 
 .SUFFIXES:.c .o
 


@@ -36,9 +35,8 @@ lib:libbam.a
 libbam.a:$(LOBJS)
                $(AR) -cru $@ $(LOBJS)
 
 libbam.a:$(LOBJS)
                $(AR) -cru $@ $(LOBJS)
 

-### For the curses library: comment out `-lcurses' if you do not have curses installed

 samtools:lib $(AOBJS)
 samtools:lib $(AOBJS)

-               $(CC) $(CFLAGS) -o $@ $(AOBJS) $(LIBPATH) -lm -lcurses -lz -L. -lbam
+               $(CC) $(CFLAGS) -o $@ $(AOBJS) -lm $(LIBPATH) $(LIBCURSES) -lz -L. -lbam

 
 razip:razip.o razf.o
                $(CC) $(CFLAGS) -o $@ razf.o razip.o -lz
 
 razip:razip.o razf.o
                $(CC) $(CFLAGS) -o $@ razf.o razip.o -lz

diff --git a/Makefile.mingw b/Makefile.mingw
new file mode 100644 (file)
index 0000000..80e8009
--- /dev/null
+++ b/Makefile.mingw
@@ -0,0 +1,55 @@
+CC=                    gcc.exe
+AR=                    ar.exe
+CFLAGS=                -g -Wall -O2
+DFLAGS=                -D_FILE_OFFSET_BITS=64 -D_CURSES_LIB=2 -D_USE_KNETFILE
+LOBJS=         bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o \
+                       bam_pileup.o bam_lpileup.o bam_md.o glf.o razf.o faidx.o bam_sort.o \
+                       knetfile.o
+AOBJS=         bam_tview.o bam_maqcns.o bam_plcmd.o sam_view.o \
+                       bam_rmdup.o bam_rmdupse.o bam_mate.o bam_stat.o bam_color.o     \
+                       bamtk.o
+PROG=          samtools
+INCLUDES=      -Iwin32
+SUBDIRS=       .
+LIBPATH=       
+
+.SUFFIXES:.c .o
+
+.c.o:
+               $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $< -o $@
+
+all:$(PROG)
+
+lib:libbam.a
+
+libbam.a:$(LOBJS)
+               $(AR) -cru $@ $(LOBJS)
+
+samtools:lib $(AOBJS)
+               $(CC) $(CFLAGS) -o $@ $(AOBJS) $(LIBPATH) -lm -L. -lbam -Lwin32 -lz -lcurses -lws2_32
+
+razip:razip.o razf.o
+               $(CC) $(CFLAGS) -o $@ razf.o razip.o -lz
+
+bgzip:bgzip.o bgzf.o
+               $(CC) $(CFLAGS) -o $@ bgzf.o bgzip.o -lz
+
+razip.o:razf.h
+bam.o:bam.h razf.h bam_endian.h kstring.h
+sam.o:sam.h bam.h
+bam_import.o:bam.h kseq.h khash.h razf.h
+bam_pileup.o:bam.h razf.h ksort.h
+bam_plcmd.o:bam.h faidx.h bam_maqcns.h glf.h
+bam_index.o:bam.h khash.h ksort.h razf.h bam_endian.h
+bam_lpileup.o:bam.h ksort.h
+bam_tview.o:bam.h faidx.h bam_maqcns.h
+bam_maqcns.o:bam.h ksort.h bam_maqcns.h
+bam_sort.o:bam.h ksort.h razf.h
+bam_md.o:bam.h faidx.h
+glf.o:glf.h
+
+faidx.o:faidx.h razf.h khash.h
+faidx_main.o:faidx.h razf.h
+
+clean:
+               rm -fr gmon.out *.o *.exe *.dSYM razip $(PROG) *~ *.a

diff --git a/NEWS b/NEWS
index 149c090c08bc9724d677dce8a66ddf62a83da66b..8e0ba351ba15c148491119e45b2c1531ded6f66c 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,54 @@
+Beta Release 0.1.6 (2 September, 2009)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Notable changes:
+
+ * In tview, do not show a blank screen when no reads mapped to the
+   corresponding region.
+
+ * Implemented native HTTP support in the BGZF library. Samtools is now
+   able to directly open a BAM file on HTTP. HTTP proxy is also
+   supported via the "http_proxy" environmental variable.
+
+ * Samtools is now compitable with the MinGW (win32) compiler and the
+   PDCurses library.
+
+ * The calmd (or fillmd) command now calculates the NM tag and replaces
+   MD tags if they are wrong.
+
+ * The view command now recognizes and optionally prints FLAG in HEXs or
+   strings to make a SAM file more friendly to human eyes. This is a
+   samtools-C extension, not implemented in Picard for the time
+   being. Please type `samtools view -?' for more information.
+
+ * BAM files now have an end-of-file (EOF) marker to facilitate
+   truncation detection. A warning will be given if an on-disk BAM file
+   does not have this marker. The warning will be seen on BAM files
+   generated by an older version of samtools. It does NO harm.
+
+ * New key bindings in tview: `r' to show read names and `s' to show
+   reference skip (N operation) as deletions.
+
+ * Fixed a bug in `samtools merge -n'.
+
+ * Samtools merge now optionally copies the header of a user specified
+   SAM file to the resultant BAM output.
+
+ * Samtools pileup/tview works with a CIGAR with the first or the last
+   operation is an indel.
+
+ * Fixed a bug in bam_aux_get().
+
+
+Changes in other utilies:
+
+ * Fixed wrong FLAG in maq2sam.
+
+
+(0.1.6: 2 September 2009, r453)
+
+
+

 Beta Release 0.1.5 (7 July, 2009)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Beta Release 0.1.5 (7 July, 2009)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

diff --git a/bam.c b/bam.c
index 1ff4a5aba25028f4fccdd0735a9dcae594b6e9a8..619b46a6718a773aa073e7ecede13ef5cd10dee8 100644 (file)
--- a/bam.c
+++ b/bam.c
@@ -6,6 +6,7 @@
 #include "kstring.h"
 
 int bam_is_be = 0;
 #include "kstring.h"
 
 int bam_is_be = 0;

+char *bam_flag2char_table = "pPuUrR12sfd\0\0\0\0\0";

 
 /**************************
  * CIGAR related routines *
 
 /**************************
  * CIGAR related routines *


@@ -90,7 +91,11 @@ bam_header_t *bam_header_read(bamFile fp)
 {
        bam_header_t *header;
        char buf[4];
 {
        bam_header_t *header;
        char buf[4];

-       int32_t i, name_len;
+       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");
+       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;
        if (strncmp(buf, "BAM\001", 4)) {
        // read "BAM1"
        if (bam_read(fp, buf, 4) != 4) return 0;
        if (strncmp(buf, "BAM\001", 4)) {


@@ -236,7 +241,7 @@ int bam_write1(bamFile fp, const bam1_t *b)
        return bam_write1_core(fp, &b->core, b->data_len, b->data);
 }
 
        return bam_write1_core(fp, &b->core, b->data_len, b->data);
 }
 

-char *bam_format1(const bam_header_t *header, const bam1_t *b)
+char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int of)

 {
        uint8_t *s = bam1_seq(b), *t = bam1_qual(b);
        int i;
 {
        uint8_t *s = bam1_seq(b), *t = bam1_qual(b);
        int i;


@@ -244,7 +249,15 @@ char *bam_format1(const bam_header_t *header, const bam1_t *b)
        kstring_t str;
        str.l = str.m = 0; str.s = 0;
 
        kstring_t str;
        str.l = str.m = 0; str.s = 0;
 

-       ksprintf(&str, "%s\t%d\t", bam1_qname(b), c->flag);
+       ksprintf(&str, "%s\t", bam1_qname(b));
+       if (of == BAM_OFDEC) ksprintf(&str, "%d\t", c->flag);
+       else if (of == BAM_OFHEX) ksprintf(&str, "0x%x\t", c->flag);
+       else { // BAM_OFSTR
+               for (i = 0; i < 16; ++i)
+                       if ((c->flag & 1<<i) && bam_flag2char_table[i])
+                               kputc(bam_flag2char_table[i], &str);
+               kputc('\t', &str);
+       }

        if (c->tid < 0) kputs("*\t", &str);
        else ksprintf(&str, "%s\t", header->target_name[c->tid]);
        ksprintf(&str, "%d\t%d\t", c->pos + 1, c->qual);
        if (c->tid < 0) kputs("*\t", &str);
        else ksprintf(&str, "%s\t", header->target_name[c->tid]);
        ksprintf(&str, "%d\t%d\t", c->pos + 1, c->qual);


@@ -258,10 +271,12 @@ char *bam_format1(const bam_header_t *header, const bam1_t *b)
        else if (c->mtid == c->tid) kputs("=\t", &str);
        else ksprintf(&str, "%s\t", header->target_name[c->mtid]);
        ksprintf(&str, "%d\t%d\t", c->mpos + 1, c->isize);
        else if (c->mtid == c->tid) kputs("=\t", &str);
        else ksprintf(&str, "%s\t", header->target_name[c->mtid]);
        ksprintf(&str, "%d\t%d\t", c->mpos + 1, c->isize);

-       for (i = 0; i < c->l_qseq; ++i) kputc(bam_nt16_rev_table[bam1_seqi(s, i)], &str);
-       kputc('\t', &str);
-       if (t[0] == 0xff) kputc('*', &str);
-       else for (i = 0; i < c->l_qseq; ++i) kputc(t[i] + 33, &str);
+       if (c->l_qseq) {
+               for (i = 0; i < c->l_qseq; ++i) kputc(bam_nt16_rev_table[bam1_seqi(s, i)], &str);
+               kputc('\t', &str);
+               if (t[0] == 0xff) kputc('*', &str);
+               else for (i = 0; i < c->l_qseq; ++i) kputc(t[i] + 33, &str);
+       } else ksprintf(&str, "*\t*");

        s = bam1_aux(b);
        while (s < b->data + b->data_len) {
                uint8_t type, key[2];
        s = bam1_aux(b);
        while (s < b->data + b->data_len) {
                uint8_t type, key[2];


@@ -282,6 +297,11 @@ char *bam_format1(const bam_header_t *header, const bam1_t *b)
        return str.s;
 }
 
        return str.s;
 }
 

+char *bam_format1(const bam_header_t *header, const bam1_t *b)
+{
+       return bam_format1_core(header, b, BAM_OFDEC);
+}
+

 void bam_view1(const bam_header_t *header, const bam1_t *b)
 {
        char *s = bam_format1(header, b);
 void bam_view1(const bam_header_t *header, const bam1_t *b)
 {
        char *s = bam_format1(header, b);

diff --git a/bam.h b/bam.h
index 83c03ada21098eda911caabb0e56db3864599e73..ec983dfdd183042ef67a019e215a6dff6724c71f 100644 (file)
--- a/bam.h
+++ b/bam.h
@@ -45,21 +45,7 @@
 #include <string.h>
 #include <stdio.h>
 
 #include <string.h>
 #include <stdio.h>
 

-#define _IOLIB 2
-
-#if _IOLIB == 1 && !defined(_NO_RAZF)
-#define BAM_TRUE_OFFSET
-#include "razf.h"
-/*! @abstract BAM file handler */
-typedef RAZF *bamFile;
-#define bam_open(fn, mode) razf_open(fn, mode)
-#define bam_dopen(fd, mode) razf_dopen(fd, mode)
-#define bam_close(fp) razf_close(fp)
-#define bam_read(fp, buf, size) razf_read(fp, buf, size)
-#define bam_write(fp, buf, size) razf_write(fp, buf, size)
-#define bam_tell(fp) razf_tell(fp)
-#define bam_seek(fp, pos, dir) razf_seek(fp, pos, dir)
-#elif _IOLIB == 2
+#ifndef BAM_LITE

 #define BAM_VIRTUAL_OFFSET16
 #include "bgzf.h"
 /*! @abstract BAM file handler */
 #define BAM_VIRTUAL_OFFSET16
 #include "bgzf.h"
 /*! @abstract BAM file handler */


@@ -71,18 +57,15 @@ typedef BGZF *bamFile;
 #define bam_write(fp, buf, size) bgzf_write(fp, buf, size)
 #define bam_tell(fp) bgzf_tell(fp)
 #define bam_seek(fp, pos, dir) bgzf_seek(fp, pos, dir)
 #define bam_write(fp, buf, size) bgzf_write(fp, buf, size)
 #define bam_tell(fp) bgzf_tell(fp)
 #define bam_seek(fp, pos, dir) bgzf_seek(fp, pos, dir)

-#elif _IOLIB == 3
-#define BAM_VIRTUAL_OFFSET16
-#include "razf.h"
-/*! @abstract BAM file handler */
-typedef RAZF *bamFile;
-#define bam_open(fn, mode) razf_open2(fn, mode)
-#define bam_dopen(fd, mode) razf_dopen2(fd, mode)
-#define bam_close(fp) razf_close(fp)
-#define bam_read(fp, buf, size) razf_read(fp, buf, size)
-#define bam_write(fp, buf, size) razf_write(fp, buf, size)
-#define bam_tell(fp) razf_tell2(fp)
-#define bam_seek(fp, pos, dir) razf_seek2(fp, pos, dir)
+#else
+#define BAM_TRUE_OFFSET
+#include <zlib.h>
+typedef gzFile bamFile;
+#define bam_open(fn, mode) gzopen(fn, mode)
+#define bam_dopen(fd, mode) gzdopen(fd, mode)
+#define bam_close(fp) gzclose(fp)
+#define bam_read(fp, buf, size) gzread(fp, buf, size)
+/* no bam_write/bam_tell/bam_seek() here */

 #endif
 
 /*! @typedef
 #endif
 
 /*! @typedef


@@ -130,6 +113,10 @@ typedef struct {
 /*! @abstract optical or PCR duplicate */
 #define BAM_FDUP        1024
 
 /*! @abstract optical or PCR duplicate */
 #define BAM_FDUP        1024
 

+#define BAM_OFDEC          0
+#define BAM_OFHEX          1
+#define BAM_OFSTR          2
+

 /*! @abstract defautl mask for pileup */
 #define BAM_DEF_MASK (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP)
 
 /*! @abstract defautl mask for pileup */
 #define BAM_DEF_MASK (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP)
 


@@ -448,6 +435,8 @@ extern "C" {
         */
        char *bam_format1(const bam_header_t *header, const bam1_t *b);
 
         */
        char *bam_format1(const bam_header_t *header, const bam1_t *b);
 

+       char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int of);
+

        /*! @typedef
          @abstract Structure for one alignment covering the pileup position.
          @field  b      pointer to the alignment
        /*! @typedef
          @abstract Structure for one alignment covering the pileup position.
          @field  b      pointer to the alignment


@@ -526,6 +515,8 @@ extern "C" {
         */
        int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf);
 
         */
        int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf);
 

+       int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data);
+

        struct __bam_lplbuf_t;
        typedef struct __bam_lplbuf_t bam_lplbuf_t;
 
        struct __bam_lplbuf_t;
        typedef struct __bam_lplbuf_t bam_lplbuf_t;
 


@@ -540,9 +531,6 @@ extern "C" {
        /*! @abstract  bam_plbuf_push() equivalent with level calculated. */
        int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *buf);
 
        /*! @abstract  bam_plbuf_push() equivalent with level calculated. */
        int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *buf);
 

-       /*! @abstract  bam_plbuf_file() equivalent with level calculated. */
-       int bam_lpileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data);
-

        struct __bam_index_t;
        typedef struct __bam_index_t bam_index_t;
 
        struct __bam_index_t;
        typedef struct __bam_index_t bam_index_t;
 


@@ -622,8 +610,8 @@ extern "C" {
        char bam_aux2A(const uint8_t *s);
        char *bam_aux2Z(const uint8_t *s);
 
        char bam_aux2A(const uint8_t *s);
        char *bam_aux2Z(const uint8_t *s);
 

+       int bam_aux_del(bam1_t *b, uint8_t *s);

        void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data);
        void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data);

-

        uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]); // an alias of bam_aux_get()
 
        /*!  
        uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]); // an alias of bam_aux_get()
 
        /*!  

diff --git a/bam_aux.c b/bam_aux.c
index 7482500ef08a53f3471858793112d9c72263e33a..d0d733fc467bfa9af7d2f032b83ffb47840ef18a 100644 (file)
--- a/bam_aux.c
+++ b/bam_aux.c
@@ -25,24 +25,41 @@ uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2])
        return bam_aux_get(b, tag);
 }
 
        return bam_aux_get(b, tag);
 }
 

+#define __skip_tag(s) do { \
+               int type = toupper(*(s));                                                                               \
+               ++(s);                                                                                                                  \
+               if (type == 'C' || type == 'A') ++(s);                                                  \
+               else if (type == 'S') (s) += 2;                                                                 \
+               else if (type == 'I' || type == 'F') (s) += 4;                                  \
+               else if (type == 'D') (s) += 8;                                                                 \
+               else if (type == 'Z' || type == 'H') { while (*(s)) ++(s); ++(s); } \
+       } while (0)
+

 uint8_t *bam_aux_get(const bam1_t *b, const char tag[2])
 {
        uint8_t *s;
        int y = tag[0]<<8 | tag[1];
        s = bam1_aux(b);
        while (s < b->data + b->data_len) {
 uint8_t *bam_aux_get(const bam1_t *b, const char tag[2])
 {
        uint8_t *s;
        int y = tag[0]<<8 | tag[1];
        s = bam1_aux(b);
        while (s < b->data + b->data_len) {

-               int type, x = (int)s[0]<<8 | s[1];
+               int x = (int)s[0]<<8 | s[1];

                s += 2;
                if (x == y) return s;
                s += 2;
                if (x == y) return s;

-               type = toupper(*s); ++s;
-               if (type == 'C') ++s;
-               else if (type == 'S') s += 2;
-               else if (type == 'I' || type == 'F') s += 4;
-               else if (type == 'D') s += 8;
-               else if (type == 'Z' || type == 'H') { while (*s) ++s; ++s; }
+               __skip_tag(s);

        }
        return 0;
 }
        }
        return 0;
 }

+// s MUST BE returned by bam_aux_get()
+int bam_aux_del(bam1_t *b, uint8_t *s)
+{
+       uint8_t *p, *aux;
+       aux = bam1_aux(b);
+       p = s - 2;
+       __skip_tag(s);
+       memmove(p, s, b->l_aux - (s - aux));
+       b->data_len -= s - p;
+       b->l_aux -= s - p;
+       return 0;
+}

 
 void bam_init_header_hash(bam_header_t *header)
 {
 
 void bam_init_header_hash(bam_header_t *header)
 {

diff --git a/bam_color.c b/bam_color.c
index 75aedd603a1798e85a40c982dd9392ee2a3e5b3f..ce637f708db56de777e86f2f0ad50478fb1a6657 100644 (file)
--- a/bam_color.c
+++ b/bam_color.c
@@ -100,8 +100,8 @@ char bam_aux_getCEi(bam1_t *b, int i)
                cs_i = strlen(cs) - 1 - i;
                // get current color
                cur_color = cs[cs_i];
                cs_i = strlen(cs) - 1 - i;
                // get current color
                cur_color = cs[cs_i];

-               // get previous base
-               prev_b = (0 == cs_i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
+               // get previous base.  Note: must rc adaptor
+               prev_b = (cs_i == 1) ? "TGCAN"[(int)bam_aux_nt2int(cs[0])] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];

                // get current base
                cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)]; 
        }
                // get current base
                cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)]; 
        }

diff --git a/bam_import.c b/bam_import.c
index fccaa022208131b27093a2b44f32e74d13a469c0..1dc906eb49f36428460be547b1582ab52432f651 100644 (file)
--- a/bam_import.c
+++ b/bam_import.c
@@ -5,6 +5,9 @@
 #include <stdlib.h>
 #include <unistd.h>
 #include <assert.h>
 #include <stdlib.h>
 #include <unistd.h>
 #include <assert.h>

+#ifdef _WIN32
+#include <fcntl.h>
+#endif

 #include "kstring.h"
 #include "bam.h"
 #include "kseq.h"
 #include "kstring.h"
 #include "bam.h"
 #include "kseq.h"


@@ -36,6 +39,25 @@ unsigned char bam_nt16_table[256] = {
        15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15
 };
 
        15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15
 };
 

+unsigned short bam_char2flag_table[256] = {
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,BAM_FREAD1,BAM_FREAD2,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       BAM_FPROPER_PAIR,0,BAM_FMREVERSE,0, 0,BAM_FMUNMAP,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, BAM_FDUP,0,BAM_FQCFAIL,0, 0,0,0,0, 0,0,0,0,
+       BAM_FPAIRED,0,BAM_FREVERSE,BAM_FSECONDARY, 0,BAM_FUNMAP,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
+       0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0
+};
+

 char *bam_nt16_rev_table = "=ACMGRSVTWYHKDBN";
 
 struct __tamFile_t {
 char *bam_nt16_rev_table = "=ACMGRSVTWYHKDBN";
 
 struct __tamFile_t {


@@ -99,9 +121,10 @@ bam_header_t *sam_header_read2(const char *fn)
        kstring_t *str;
        kh_ref_t *hash;
        khiter_t k;
        kstring_t *str;
        kh_ref_t *hash;
        khiter_t k;

-       hash = kh_init(ref);
+       if (fn == 0) return 0;

        fp = (strcmp(fn, "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(fn, "r");
        fp = (strcmp(fn, "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(fn, "r");

-       assert(fp);
+       if (fp == 0) return 0;
+       hash = kh_init(ref);

        ks = ks_init(fp);
        str = (kstring_t*)calloc(1, sizeof(kstring_t));
        while (ks_getuntil(ks, 0, str, &dret) > 0) {
        ks = ks_init(fp);
        str = (kstring_t*)calloc(1, sizeof(kstring_t));
        while (ks_getuntil(ks, 0, str, &dret) > 0) {


@@ -296,8 +319,19 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b)
                memcpy(alloc_data(b, doff + c->l_qname) + doff, str->s, c->l_qname);
                doff += c->l_qname;
        }
                memcpy(alloc_data(b, doff + c->l_qname) + doff, str->s, c->l_qname);
                doff += c->l_qname;
        }

-       { // flag, tid, pos, qual
-               ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; c->flag = atoi(str->s);
+       { // flag
+               long flag;
+               char *s;
+               ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1;
+               flag = strtol((char*)str->s, &s, 0);
+               if (*s) { // not the end of the string
+                       flag = 0;
+                       for (s = str->s; *s; ++s)
+                               flag |= bam_char2flag_table[(int)*s];
+               }
+               c->flag = flag;
+       }
+       { // tid, pos, qual

                ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; c->tid = bam_get_tid(header, str->s);
                if (c->tid < 0 && strcmp(str->s, "*")) {
                        if (header->n_targets == 0) {
                ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; c->tid = bam_get_tid(header, str->s);
                if (c->tid < 0 && strcmp(str->s, "*")) {
                        if (header->n_targets == 0) {


@@ -325,7 +359,7 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b)
                        for (i = 0, s = str->s; i != c->n_cigar; ++i) {
                                x = strtol(s, &t, 10);
                                op = toupper(*t);
                        for (i = 0, s = str->s; i != c->n_cigar; ++i) {
                                x = strtol(s, &t, 10);
                                op = toupper(*t);

-                               if (op == 'M') op = BAM_CMATCH;
+                               if (op == 'M' || op == '=' || op == 'X') op = BAM_CMATCH;

                                else if (op == 'I') op = BAM_CINS;
                                else if (op == 'D') op = BAM_CDEL;
                                else if (op == 'N') op = BAM_CREF_SKIP;
                                else if (op == 'I') op = BAM_CINS;
                                else if (op == 'D') op = BAM_CDEL;
                                else if (op == 'N') op = BAM_CREF_SKIP;


@@ -339,7 +373,13 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b)
                        if (*s) parse_error(fp->n_lines, "unmatched CIGAR operation");
                        c->bin = bam_reg2bin(c->pos, bam_calend(c, bam1_cigar(b)));
                        doff += c->n_cigar * 4;
                        if (*s) parse_error(fp->n_lines, "unmatched CIGAR operation");
                        c->bin = bam_reg2bin(c->pos, bam_calend(c, bam1_cigar(b)));
                        doff += c->n_cigar * 4;

-               } else c->bin = bam_reg2bin(c->pos, c->pos + 1);
+               } else {
+                       if (!(c->flag&BAM_FUNMAP)) {
+                               fprintf(stderr, "Parse warning at line %lld: mapped sequence without CIGAR\n", (long long)fp->n_lines);
+                               c->flag |= BAM_FUNMAP;
+                       }
+                       c->bin = bam_reg2bin(c->pos, c->pos + 1);
+               }

        }
        { // mtid, mpos, isize
                ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1;
        }
        { // mtid, mpos, isize
                ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1;


@@ -352,16 +392,18 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b)
        }
        { // seq and qual
                int i;
        }
        { // seq and qual
                int i;

-               uint8_t *p;
+               uint8_t *p = 0;

                if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -5; // seq
                z += str->l + 1;
                if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -5; // seq
                z += str->l + 1;

-               c->l_qseq = strlen(str->s);
-               if (c->n_cigar && c->l_qseq != (int32_t)bam_cigar2qlen(c, bam1_cigar(b)))
-                       parse_error(fp->n_lines, "CIGAR and sequence length are inconsistent");
-               p = (uint8_t*)alloc_data(b, doff + c->l_qseq + (c->l_qseq+1)/2) + doff;
-               bzero(p, (c->l_qseq+1)/2);
-               for (i = 0; i < c->l_qseq; ++i)
-                       p[i/2] |= bam_nt16_table[(int)str->s[i]] << 4*(1-i%2);
+               if (strcmp(str->s, "*")) {
+                       c->l_qseq = strlen(str->s);
+                       if (c->n_cigar && c->l_qseq != (int32_t)bam_cigar2qlen(c, bam1_cigar(b)))
+                               parse_error(fp->n_lines, "CIGAR and sequence length are inconsistent");
+                       p = (uint8_t*)alloc_data(b, doff + c->l_qseq + (c->l_qseq+1)/2) + doff;
+                       memset(p, 0, (c->l_qseq+1)/2);
+                       for (i = 0; i < c->l_qseq; ++i)
+                               p[i/2] |= bam_nt16_table[(int)str->s[i]] << 4*(1-i%2);
+               } else c->l_qseq = 0;

                if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -6; // qual
                z += str->l + 1;
                if (strcmp(str->s, "*") && c->l_qseq != strlen(str->s))
                if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -6; // qual
                z += str->l + 1;
                if (strcmp(str->s, "*") && c->l_qseq != strlen(str->s))


@@ -457,9 +499,11 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b)
 tamFile sam_open(const char *fn)
 {
        tamFile fp;
 tamFile sam_open(const char *fn)
 {
        tamFile fp;

+       gzFile gzfp = (strcmp(fn, "-") == 0)? gzdopen(fileno(stdin), "rb") : gzopen(fn, "rb");
+       if (gzfp == 0) return 0;

        fp = (tamFile)calloc(1, sizeof(struct __tamFile_t));
        fp->str = (kstring_t*)calloc(1, sizeof(kstring_t));
        fp = (tamFile)calloc(1, sizeof(struct __tamFile_t));
        fp->str = (kstring_t*)calloc(1, sizeof(kstring_t));

-       fp->fp = (strcmp(fn, "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(fn, "r");
+       fp->fp = gzfp;

        fp->ks = ks_init(fp->fp);
        return fp;
 }
        fp->ks = ks_init(fp->fp);
        return fp;
 }

diff --git a/bam_index.c b/bam_index.c
index 72ef270f9c7afd143759e1f98d681ce69147bb13..4ff6bd482505172741c5a48ab413f141342e93bc 100644 (file)
--- a/bam_index.c
+++ b/bam_index.c
@@ -4,7 +4,9 @@
 #include "khash.h"
 #include "ksort.h"
 #include "bam_endian.h"
 #include "khash.h"
 #include "ksort.h"
 #include "bam_endian.h"

+#ifdef _USE_KNETFILE

 #include "knetfile.h"
 #include "knetfile.h"

+#endif

 
 /*!
   @header
 
 /*!
   @header


@@ -328,7 +330,7 @@ bam_index_t *bam_index_load_local(const char *_fn)
        FILE *fp;
        char *fnidx, *fn;
 
        FILE *fp;
        char *fnidx, *fn;
 

-       if (strstr(_fn, "ftp://") == _fn) {
+       if (strstr(_fn, "ftp://") == _fn || strstr(_fn, "http://") == _fn) {

                const char *p;
                int l = strlen(_fn);
                for (p = _fn + l - 1; p >= _fn; --p)
                const char *p;
                int l = strlen(_fn);
                for (p = _fn + l - 1; p >= _fn; --p)


@@ -354,6 +356,7 @@ bam_index_t *bam_index_load_local(const char *_fn)
        } else return 0;
 }
 
        } else return 0;
 }
 

+#ifdef _USE_KNETFILE

 static void download_from_remote(const char *url)
 {
        const int buf_size = 1 * 1024 * 1024;
 static void download_from_remote(const char *url)
 {
        const int buf_size = 1 * 1024 * 1024;


@@ -362,7 +365,7 @@ static void download_from_remote(const char *url)
        uint8_t *buf;
        knetFile *fp_remote;
        int l;
        uint8_t *buf;
        knetFile *fp_remote;
        int l;

-       if (strstr(url, "ftp://") != url) return;
+       if (strstr(url, "ftp://") != url && strstr(url, "http://") != url) return;

        l = strlen(url);
        for (fn = (char*)url + l - 1; fn >= url; --fn)
                if (*fn == '/') break;
        l = strlen(url);
        for (fn = (char*)url + l - 1; fn >= url; --fn)
                if (*fn == '/') break;


@@ -384,12 +387,18 @@ static void download_from_remote(const char *url)
        fclose(fp);
        knet_close(fp_remote);
 }
        fclose(fp);
        knet_close(fp_remote);
 }

+#else
+static void download_from_remote(const char *url)
+{
+       return;
+}
+#endif

 
 bam_index_t *bam_index_load(const char *fn)
 {
        bam_index_t *idx;
        idx = bam_index_load_local(fn);
 
 bam_index_t *bam_index_load(const char *fn)
 {
        bam_index_t *idx;
        idx = bam_index_load_local(fn);

-       if (idx == 0 && strstr(fn, "ftp://") == fn) {
+       if (idx == 0 && (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn)) {

                char *fnidx = calloc(strlen(fn) + 5, 1);
                strcat(strcpy(fnidx, fn), ".bai");
                fprintf(stderr, "[bam_index_load] attempting to download the remote index file.\n");
                char *fnidx = calloc(strlen(fn) + 5, 1);
                strcat(strcpy(fnidx, fn), ".bai");
                fprintf(stderr, "[bam_index_load] attempting to download the remote index file.\n");

diff --git a/bam_lpileup.c b/bam_lpileup.c
index 425290e26f11d29cd9b3cc4633287ca5e8ea41e2..d4dd63b99b85dd4c5811758153a56757c908a52c 100644 (file)
--- a/bam_lpileup.c
+++ b/bam_lpileup.c
@@ -196,19 +196,3 @@ int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *tv)
 {
        return bam_plbuf_push(b, tv->plbuf);
 }
 {
        return bam_plbuf_push(b, tv->plbuf);
 }

-
-int bam_lpileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data)
-{
-       bam_lplbuf_t *buf;
-       int ret;
-       bam1_t *b;
-       b = (bam1_t*)calloc(1, sizeof(bam1_t));
-       buf = bam_lplbuf_init(func, func_data);
-       bam_plbuf_set_mask(buf->plbuf, mask);
-       while ((ret = bam_read1(fp, b)) >= 0)
-               bam_lplbuf_push(b, buf);
-       bam_lplbuf_push(0, buf);
-       bam_lplbuf_destroy(buf);
-       free(b->data); free(b);
-       return 0;
-}

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 464288ae7da0ac74f1443bbb55c03684927cb391..f36b0ee2ab443affe0635866a8d593c5cb54fdf7 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -439,7 +439,7 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                }
                                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 (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;
+                               //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]);
                        }
                                //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]);
                        }


@@ -499,6 +499,15 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                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;
                        }
                                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 cnt_ref and cnt_ambi
+                       if (max1_i != 0 && max2_i != 0) {
+                               for (j = 0; j < n; ++j) {
+                                       int diff1 = score[j] - score[max1_i * n + j];
+                                       int diff2 = score[j] - score[max2_i * n + j];
+                                       if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
+                                       else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
+                               }
+                       }

                }
                free(score); free(pscore); free(ref2); free(inscns);
        }
                }
                free(score); free(pscore); free(ref2); free(inscns);
        }

diff --git a/bam_maqcns.h b/bam_maqcns.h
index 36704d7bce51104e1a8d73ed233f8ebd4b8ce582..2a82aeee8057f5992570d0896a1621b9116057cd 100644 (file)
--- a/bam_maqcns.h
+++ b/bam_maqcns.h
@@ -24,7 +24,8 @@ typedef struct {
 
 typedef struct {
        int indel1, indel2;
 
 typedef struct {
        int indel1, indel2;

-       int cnt1, cnt2, cnt_ambi, cnt_anti;
+       int cnt1, cnt2, cnt_anti;
+       int cnt_ref, cnt_ambi;

        char *s[2];
        //
        int gt, gl[2];
        char *s[2];
        //
        int gt, gl[2];

diff --git a/bam_md.c b/bam_md.c
index a20f9b354d477c1f11a8d98854ff9ec7682d01c1..8d074870d336677cee2e606ab05ee81a2b4a6f1e 100644 (file)
--- a/bam_md.c
+++ b/bam_md.c
@@ -3,7 +3,7 @@
 #include <string.h>
 #include <ctype.h>
 #include "faidx.h"
 #include <string.h>
 #include <ctype.h>
 #include "faidx.h"

-#include "bam.h"
+#include "sam.h"

 #include "kstring.h"
 
 void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
 #include "kstring.h"
 
 void bam_fillmd1(bam1_t *b, char *ref, int is_equal)


@@ -13,11 +13,9 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
        bam1_core_t *c = &b->core;
        int i, x, y, u = 0;
        kstring_t *str;
        bam1_core_t *c = &b->core;
        int i, x, y, u = 0;
        kstring_t *str;

-       uint8_t *old_md;
+       uint8_t *old_md, *old_nm;
+       int32_t old_nm_i = -1, nm = 0;

 
 

-       old_md = bam_aux_get(b, "MD");
-       if (c->flag & BAM_FUNMAP) return;
-       if (old_md && !is_equal) return; // no need to add MD

        str = (kstring_t*)calloc(1, sizeof(kstring_t));
        for (i = y = 0, x = c->pos; i < c->n_cigar; ++i) {
                int j, l = cigar[i]>>4, op = cigar[i]&0xf;
        str = (kstring_t*)calloc(1, sizeof(kstring_t));
        for (i = y = 0, x = c->pos; i < c->n_cigar; ++i) {
                int j, l = cigar[i]>>4, op = cigar[i]&0xf;


@@ -26,13 +24,13 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
                                int z = y + j;
                                int c1 = bam1_seqi(seq, z), c2 = bam_nt16_table[(int)ref[x+j]];
                                if (ref[x+j] == 0) break; // out of boundary
                                int z = y + j;
                                int c1 = bam1_seqi(seq, z), c2 = bam_nt16_table[(int)ref[x+j]];
                                if (ref[x+j] == 0) break; // out of boundary

-                               if ((c1 == c2 && c1 != 15 && c2 != 15) || c1 == 0) {
+                               if ((c1 == c2 && c1 != 15 && c2 != 15) || c1 == 0) { // a match

                                        if (is_equal) seq[z/2] &= (z&1)? 0xf0 : 0x0f;
                                        ++u;
                                } else {
                                        ksprintf(str, "%d", u);
                                        kputc(ref[x+j], str);
                                        if (is_equal) seq[z/2] &= (z&1)? 0xf0 : 0x0f;
                                        ++u;
                                } else {
                                        ksprintf(str, "%d", u);
                                        kputc(ref[x+j], str);

-                                       u = 0;
+                                       u = 0; ++nm;

                                }
                        }
                        if (j < l) break;
                                }
                        }
                        if (j < l) break;


@@ -46,14 +44,27 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
                        }
                        u = 0;
                        if (j < l) break;
                        }
                        u = 0;
                        if (j < l) break;

-                       x += l;
+                       x += l; nm += l;

                } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
                        y += l;
                } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
                        y += l;

+                       if (op == BAM_CINS) nm += l;

                } else if (op == BAM_CREF_SKIP) {
                        x += l;
                }
        }
        ksprintf(str, "%d", u);
                } else if (op == BAM_CREF_SKIP) {
                        x += l;
                }
        }
        ksprintf(str, "%d", u);

+       // update NM
+       old_nm = bam_aux_get(b, "NM");
+       if (c->flag & BAM_FUNMAP) return;
+       if (old_nm) old_nm_i = bam_aux2i(old_nm);
+       if (!old_nm) bam_aux_append(b, "NM", 'i', 4, (uint8_t*)&nm);
+       else if (nm != old_nm_i) {
+               fprintf(stderr, "[bam_fillmd1] different NM for read '%s': %d -> %d\n", bam1_qname(b), old_nm_i, nm);
+               bam_aux_del(b, old_nm);
+               bam_aux_append(b, "NM", 'i', 4, (uint8_t*)&nm);
+       }
+       // update MD
+       old_md = bam_aux_get(b, "MD");

        if (!old_md) bam_aux_append(b, "MD", 'Z', str->l + 1, (uint8_t*)str->s);
        else {
                int is_diff = 0;
        if (!old_md) bam_aux_append(b, "MD", 'Z', str->l + 1, (uint8_t*)str->s);
        else {
                int is_diff = 0;


@@ -63,55 +74,73 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
                                        break;
                        if (i < str->l) is_diff = 1;
                } else is_diff = 1;
                                        break;
                        if (i < str->l) is_diff = 1;
                } else is_diff = 1;

-               if (is_diff)
-                       fprintf(stderr, "[bam_fillmd1] different MD for read '%s': '%s' != '%s'\n", bam1_qname(b), old_md+1, str->s);
+               if (is_diff) {
+                       fprintf(stderr, "[bam_fillmd1] different MD for read '%s': '%s' -> '%s'\n", bam1_qname(b), old_md+1, str->s);
+                       bam_aux_del(b, old_md);
+                       bam_aux_append(b, "MD", 'Z', str->l + 1, (uint8_t*)str->s);
+               }

        }
        free(str->s); free(str);
 }
 
 int bam_fillmd(int argc, char *argv[])
 {
        }
        free(str->s); free(str);
 }
 
 int bam_fillmd(int argc, char *argv[])
 {

-       int c, is_equal = 0, tid = -2, ret, len;
-       bamFile fp, fpout = 0;
-       bam_header_t *header;
+       int c, is_equal = 0, tid = -2, ret, len, is_bam_out, is_sam_in, is_uncompressed;
+       samfile_t *fp, *fpout = 0;

        faidx_t *fai;
        faidx_t *fai;

-       char *ref = 0;
+       char *ref = 0, mode_w[8], mode_r[8];

        bam1_t *b;
 
        bam1_t *b;
 

-       while ((c = getopt(argc, argv, "e")) >= 0) {
+       is_bam_out = is_sam_in = is_uncompressed = 0;
+       mode_w[0] = mode_r[0] = 0;
+       strcpy(mode_r, "r"); strcpy(mode_w, "w");
+       while ((c = getopt(argc, argv, "eubS")) >= 0) {

                switch (c) {
                case 'e': is_equal = 1; break;
                switch (c) {
                case 'e': is_equal = 1; break;

+               case 'b': is_bam_out = 1; break;
+               case 'u': is_uncompressed = is_bam_out = 1; break;
+               case 'S': is_sam_in = 1; break;

                default: fprintf(stderr, "[bam_fillmd] unrecognized option '-%c'\n", c); return 1;
                }
        }
                default: fprintf(stderr, "[bam_fillmd] unrecognized option '-%c'\n", c); return 1;
                }
        }

+       if (!is_sam_in) strcat(mode_r, "b");
+       if (is_bam_out) strcat(mode_w, "b");
+       else strcat(mode_w, "h");
+       if (is_uncompressed) strcat(mode_w, "u");

        if (optind + 1 >= argc) {
        if (optind + 1 >= argc) {

-               fprintf(stderr, "Usage: bam fillmd [-e] <aln.bam> <ref.fasta>\n");
+               fprintf(stderr, "\n");
+               fprintf(stderr, "Usage:   samtools fillmd [-eubS] <aln.bam> <ref.fasta>\n\n");
+               fprintf(stderr, "Options: -e       change identical bases to '='\n");
+               fprintf(stderr, "         -u       uncompressed BAM output (for piping)\n");
+               fprintf(stderr, "         -b       compressed BAM output\n");
+               fprintf(stderr, "         -S       the input is SAM with header\n\n");
+               return 1;
+       }
+       fp = samopen(argv[optind], mode_r, 0);
+       if (fp == 0) return 1;
+       if (is_sam_in && (fp->header == 0 || fp->header->n_targets == 0)) {
+               fprintf(stderr, "[bam_fillmd] input SAM does not have header. Abort!\n");

                return 1;
        }
                return 1;
        }

-       fp = strcmp(argv[optind], "-")? bam_open(argv[optind], "r") : bam_dopen(fileno(stdin), "r");
-       assert(fp);
-       header = bam_header_read(fp);
-       fpout = bam_dopen(fileno(stdout), "w");
-       bam_header_write(fpout, header);
+       fpout = samopen("-", mode_w, fp->header);

        fai = fai_load(argv[optind+1]);
 
        b = bam_init1();
        fai = fai_load(argv[optind+1]);
 
        b = bam_init1();

-       while ((ret = bam_read1(fp, b)) >= 0) {
+       while ((ret = samread(fp, b)) >= 0) {

                if (b->core.tid >= 0) {
                        if (tid != b->core.tid) {
                                free(ref);
                if (b->core.tid >= 0) {
                        if (tid != b->core.tid) {
                                free(ref);

-                               ref = fai_fetch(fai, header->target_name[b->core.tid], &len);
+                               ref = fai_fetch(fai, fp->header->target_name[b->core.tid], &len);

                                tid = b->core.tid;
                        }
                        bam_fillmd1(b, ref, is_equal);
                }
                                tid = b->core.tid;
                        }
                        bam_fillmd1(b, ref, is_equal);
                }

-               bam_write1(fpout, b);
+               samwrite(fpout, b);

        }
        bam_destroy1(b);
 
        free(ref);
        fai_destroy(fai);
        }
        bam_destroy1(b);
 
        free(ref);
        fai_destroy(fai);

-       bam_header_destroy(header);
-       bam_close(fp); bam_close(fpout);
+       samclose(fp); samclose(fpout);

        return 0;
 }
        return 0;
 }

diff --git a/bam_pileup.c b/bam_pileup.c
index 3ffd52851e3be64bed628081b2b57d8ca6378b8b..f68f400ab5cd086dc88a33e14899abe6cd4e53bb 100644 (file)
--- a/bam_pileup.c
+++ b/bam_pileup.c
@@ -78,6 +78,10 @@ static inline int resolve_cigar(bam_pileup1_t *p, uint32_t pos)
                                                int op_next = cigar&BAM_CIGAR_MASK; // next CIGAR operation
                                                if (op_next == BAM_CDEL) p->indel = -(int32_t)(cigar>>BAM_CIGAR_SHIFT); // del
                                                else if (op_next == BAM_CINS) p->indel = cigar>>BAM_CIGAR_SHIFT; // ins
                                                int op_next = cigar&BAM_CIGAR_MASK; // next CIGAR operation
                                                if (op_next == BAM_CDEL) p->indel = -(int32_t)(cigar>>BAM_CIGAR_SHIFT); // del
                                                else if (op_next == BAM_CINS) p->indel = cigar>>BAM_CIGAR_SHIFT; // ins

+                                               if (op_next == BAM_CDEL || op_next == BAM_CINS) {
+                                                       if (k + 2 < c->n_cigar) op_next = bam1_cigar(b)[k+2]&BAM_CIGAR_MASK;
+                                                       else p->is_tail = 1;
+                                               }

                                                if (op_next == BAM_CSOFT_CLIP || op_next == BAM_CREF_SKIP || op_next == BAM_CHARD_CLIP)
                                                        p->is_tail = 1; // tail
                                        } else p->is_tail = 1; // this is the last operation; set tail
                                                if (op_next == BAM_CSOFT_CLIP || op_next == BAM_CREF_SKIP || op_next == BAM_CHARD_CLIP)
                                                        p->is_tail = 1; // tail
                                        } else p->is_tail = 1; // this is the last operation; set tail


@@ -166,9 +170,13 @@ int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf)
                if (b->core.flag & buf->flag_mask) return 0;
                bam_copy1(&buf->tail->b, b);
                buf->tail->beg = b->core.pos; buf->tail->end = bam_calend(&b->core, bam1_cigar(b));
                if (b->core.flag & buf->flag_mask) return 0;
                bam_copy1(&buf->tail->b, b);
                buf->tail->beg = b->core.pos; buf->tail->end = bam_calend(&b->core, bam1_cigar(b));

-               if (!(b->core.tid >= buf->max_tid || (b->core.tid == buf->max_tid && buf->tail->beg >= buf->max_pos))) {
-                       fprintf(stderr, "[bam_pileup_core] the input is not sorted. Abort!\n");
-                       abort();
+               if (b->core.tid < buf->max_tid) {
+                       fprintf(stderr, "[bam_pileup_core] the input is not sorted (chromosomes out of order)\n");
+                       return -1;
+               }
+               if ((b->core.tid == buf->max_tid) && (buf->tail->beg < buf->max_pos)) {
+                       fprintf(stderr, "[bam_pileup_core] the input is not sorted (reads out of order)\n");
+                       return -1;

                }
                buf->max_tid = b->core.tid; buf->max_pos = buf->tail->beg;
                if (buf->tail->end > buf->pos || buf->tail->b.core.tid > buf->tid) {
                }
                buf->max_tid = b->core.tid; buf->max_pos = buf->tail->beg;
                if (buf->tail->end > buf->pos || buf->tail->b.core.tid > buf->tid) {


@@ -212,3 +220,19 @@ int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf)
        }
        return 0;
 }
        }
        return 0;
 }

+
+int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data)
+{
+       bam_plbuf_t *buf;
+       int ret;
+       bam1_t *b;
+       b = bam_init1();
+       buf = bam_plbuf_init(func, func_data);
+       bam_plbuf_set_mask(buf, mask);
+       while ((ret = bam_read1(fp, b)) >= 0)
+               bam_plbuf_push(b, buf);
+       bam_plbuf_push(0, buf);
+       bam_plbuf_destroy(buf);
+       bam_destroy1(b);
+       return 0;
+}

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 5d5506fb5630a088f692e6f0c64e8b827f6a5329..5bf1ed095096b7cd91fc99d2256467b89e394d33 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -284,7 +284,8 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p
                printf("%d\t%d\t", rms_mapq, n);
                printf("%s\t%s\t", r->s[0], r->s[1]);
                //printf("%d\t%d\t", r->gl[0], r->gl[1]);
                printf("%d\t%d\t", rms_mapq, n);
                printf("%s\t%s\t", r->s[0], r->s[1]);
                //printf("%d\t%d\t", r->gl[0], r->gl[1]);

-               printf("%d\t%d\t%d\n", r->cnt1, r->cnt2, r->cnt_anti);
+               printf("%d\t%d\t%d\t", r->cnt1, r->cnt2, r->cnt_anti);
+               printf("%d\t%d\n", r->cnt_ref, r->cnt_ambi);

                bam_maqindel_ret_destroy(r);
        }
        return 0;
                bam_maqindel_ret_destroy(r);
        }
        return 0;


@@ -330,7 +331,7 @@ int bam_pileup(int argc, char *argv[])
                fprintf(stderr, "        -i        only show lines/consensus with indels\n");
                fprintf(stderr, "        -m INT    filtering reads with bits in INT [%d]\n", d->mask);
                fprintf(stderr, "        -M INT    cap mapping quality at INT [%d]\n", d->c->cap_mapQ);
                fprintf(stderr, "        -i        only show lines/consensus with indels\n");
                fprintf(stderr, "        -m INT    filtering reads with bits in INT [%d]\n", d->mask);
                fprintf(stderr, "        -M INT    cap mapping quality at INT [%d]\n", d->c->cap_mapQ);

-               fprintf(stderr, "        -t FILE   list of reference sequences (assume the input is in SAM)\n");
+               fprintf(stderr, "        -t FILE   list of reference sequences (force -S)\n");

                fprintf(stderr, "        -l FILE   list of sites at which pileup is output\n");
                fprintf(stderr, "        -f FILE   reference sequence in the FASTA format\n\n");
                fprintf(stderr, "        -c        output the maq consensus sequence\n");
                fprintf(stderr, "        -l FILE   list of sites at which pileup is output\n");
                fprintf(stderr, "        -f FILE   reference sequence in the FASTA format\n\n");
                fprintf(stderr, "        -c        output the maq consensus sequence\n");


@@ -356,6 +357,8 @@ int bam_pileup(int argc, char *argv[])
        }
        if (d->fai == 0 && (d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)))
                fprintf(stderr, "[bam_pileup] indels will not be called when -f is absent.\n");
        }
        if (d->fai == 0 && (d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)))
                fprintf(stderr, "[bam_pileup] indels will not be called when -f is absent.\n");

+       if (fn_fa && is_SAM && fn_list == 0) fn_list = samfaipath(fn_fa);
+

        {
                samfile_t *fp;
                fp = is_SAM? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0);
        {
                samfile_t *fp;
                fp = is_SAM? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0);

diff --git a/bam_rmdup.c b/bam_rmdup.c
index 1fa6cad5d7582b9b7a1d148c47df4e00aa5b638e..5da946011f10ec72758c18cc10af75180a7d8f4e 100644 (file)
--- a/bam_rmdup.c
+++ b/bam_rmdup.c
@@ -128,7 +128,8 @@ int bam_rmdup(int argc, char *argv[])
 {
        bamFile in, out;
        if (argc < 3) {
 {
        bamFile in, out;
        if (argc < 3) {

-               fprintf(stderr, "Usage: samtools rmdup <input.srt.bam> <output.bam>\n");
+               fprintf(stderr, "Usage: samtools rmdup <input.srt.bam> <output.bam>\n\n");
+               fprintf(stderr, "Note: Picard is recommended for this task.\n");

                return 1;
        }
        in = (strcmp(argv[1], "-") == 0)? bam_dopen(fileno(stdin), "r") : bam_open(argv[1], "r");
                return 1;
        }
        in = (strcmp(argv[1], "-") == 0)? bam_dopen(fileno(stdin), "r") : bam_open(argv[1], "r");

diff --git a/bam_rmdupse.c b/bam_rmdupse.c
index df03717232b10f9e23abe8ae0e7eefbc5a6c7d02..cf1b7bd7b71461938bf43033b8c452eb4a914eaf 100644 (file)
--- a/bam_rmdupse.c
+++ b/bam_rmdupse.c
@@ -161,7 +161,8 @@ int bam_rmdupse(int argc, char *argv[])
        samfile_t *in, *out;
        buffer_t *buf;
        if (argc < 3) {
        samfile_t *in, *out;
        buffer_t *buf;
        if (argc < 3) {

-               fprintf(stderr, "Usage: samtools rmdupse <in.bam> <out.bam>\n");
+               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));
                return 1;
        }
        buf = calloc(1, sizeof(buffer_t));

diff --git a/bam_sort.c b/bam_sort.c
index 402792aea44952a5962e340e9a0efbad75654494..a2d3d09473730715b78f38d6abb2b0123589daca 100644 (file)
--- a/bam_sort.c
+++ b/bam_sort.c
@@ -40,30 +40,53 @@ typedef struct {
 static inline int heap_lt(const heap1_t a, const heap1_t b)
 {
        if (g_is_by_qname) {
 static inline int heap_lt(const heap1_t a, const heap1_t b)
 {
        if (g_is_by_qname) {

-               int t = strnum_cmp(bam1_qname(a.b), bam1_qname(b.b));
+               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);
 }
 
 KSORT_INIT(heap, heap1_t, heap_lt)
 
                return (t > 0 || (t == 0 && a.pos > b.pos));
        } else return (a.pos > b.pos);
 }
 
 KSORT_INIT(heap, heap1_t, heap_lt)
 

+static void swap_header_text(bam_header_t *h1, bam_header_t *h2)
+{
+       int tempi;
+       char *temps;
+       tempi = h1->l_text, h1->l_text = h2->l_text, h2->l_text = tempi;
+       temps = h1->text, h1->text = h2->text, h2->text = temps;
+}
+

 /*!
   @abstract    Merge multiple sorted BAM.
   @param  is_by_qname whether to sort by query name
   @param  out  output BAM file name
 /*!
   @abstract    Merge multiple sorted BAM.
   @param  is_by_qname whether to sort by query name
   @param  out  output BAM file name

+  @param  headers  name of SAM file from which to copy '@' header lines,
+                   or NULL to copy them from the first file to be merged

   @param  n    number of files to be merged
   @param  fn   names of files to be merged
 
   @discussion Padding information may NOT correctly maintained. This
   function is NOT thread safe.
  */
   @param  n    number of files to be merged
   @param  fn   names of files to be merged
 
   @discussion Padding information may NOT correctly maintained. This
   function is NOT thread safe.
  */

-void bam_merge_core(int by_qname, const char *out, int n, char * const *fn)
+void bam_merge_core(int by_qname, const char *out, const char *headers, int n, char * const *fn)

 {
        bamFile fpout, *fp;
        heap1_t *heap;
        bam_header_t *hout = 0;
 {
        bamFile fpout, *fp;
        heap1_t *heap;
        bam_header_t *hout = 0;

+       bam_header_t *hheaders = NULL;

        int i, j;
 
        int i, j;
 

+       if (headers) {
+               tamFile fpheaders = sam_open(headers);
+               if (fpheaders == 0) {
+                       fprintf(stderr, "[bam_merge_core] Cannot open file `%s'. Continue anyway.\n", headers);
+               } else {
+                       hheaders = sam_header_read(fpheaders);
+                       sam_close(fpheaders);
+               }
+       }
+

        g_is_by_qname = by_qname;
        fp = (bamFile*)calloc(n, sizeof(bamFile));
        heap = (heap1_t*)calloc(n, sizeof(heap1_t));
        g_is_by_qname = by_qname;
        fp = (bamFile*)calloc(n, sizeof(bamFile));
        heap = (heap1_t*)calloc(n, sizeof(heap1_t));


@@ -72,8 +95,24 @@ void bam_merge_core(int by_qname, const char *out, int n, char * const *fn)
                bam_header_t *hin;
                assert(fp[i] = bam_open(fn[i], "r"));
                hin = bam_header_read(fp[i]);
                bam_header_t *hin;
                assert(fp[i] = bam_open(fn[i], "r"));
                hin = bam_header_read(fp[i]);

-               if (i == 0) hout = hin;
-               else { // validate multiple baf
+               if (i == 0) { // the first SAM
+                       hout = hin;
+                       if (hheaders) {
+                               // If the text headers to be swapped in include any @SQ headers,
+                               // check that they are consistent with the existing binary list
+                               // of reference information.
+                               if (hheaders->n_targets > 0) {
+                                       if (hout->n_targets != hheaders->n_targets)
+                                               fprintf(stderr, "[bam_merge_core] number of @SQ headers in `%s' differs from number of target sequences", headers);
+                                       for (j = 0; j < hout->n_targets; ++j)
+                                               if (strcmp(hout->target_name[j], hheaders->target_name[j]) != 0)
+                                                       fprintf(stderr, "[bam_merge_core] @SQ header '%s' in '%s' differs from target sequence", hheaders->target_name[j], headers);
+                               }
+                               swap_header_text(hout, hheaders);
+                               bam_header_destroy(hheaders);
+                               hheaders = NULL;
+                       }
+               } else { // validate multiple baf

                        if (hout->n_targets != hin->n_targets) {
                                fprintf(stderr, "[bam_merge_core] file '%s' has different number of target sequences. Abort!\n", fn[i]);
                                exit(1);
                        if (hout->n_targets != hin->n_targets) {
                                fprintf(stderr, "[bam_merge_core] file '%s' has different number of target sequences. Abort!\n", fn[i]);
                                exit(1);


@@ -84,9 +123,6 @@ void bam_merge_core(int by_qname, const char *out, int n, char * const *fn)
                                                        hout->target_name[j], hin->target_name[j], fn[i]);
                                        exit(1);
                                }
                                                        hout->target_name[j], hin->target_name[j], fn[i]);
                                        exit(1);
                                }

-                               if (hout->target_len[j] != hin->target_len[j])
-                                       fprintf(stderr, "[bam_merge_core] different target sequence length: %d != %d in file '%s'. Continue.\n",
-                                                       hout->target_len[j], hin->target_len[j], fn[i]);

                        }
                        bam_header_destroy(hin);
                }
                        }
                        bam_header_destroy(hin);
                }


@@ -106,34 +142,43 @@ void bam_merge_core(int by_qname, const char *out, int n, char * const *fn)
        while (heap->pos != HEAP_EMPTY) {
                bam1_t *b = heap->b;
                bam_write1_core(fpout, &b->core, b->data_len, b->data);
        while (heap->pos != HEAP_EMPTY) {
                bam1_t *b = heap->b;
                bam_write1_core(fpout, &b->core, b->data_len, b->data);

-               if ((j = bam_read1(fp[heap->i], b)) >= 0)
+               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->pos = ((uint64_t)b->core.tid<<32) | (uint32_t)b->core.pos<<1 | bam1_strand(b);

-               else if (j == -1) heap->pos = HEAP_EMPTY;
-               else fprintf(stderr, "[bam_merge_core] '%s' is truncated. Continue anyway.\n", fn[heap->i]);
+               } else if (j == -1) {
+                       heap->pos = HEAP_EMPTY;
+                       free(heap->b->data); free(heap->b);
+                       heap->b = 0;
+               } else fprintf(stderr, "[bam_merge_core] '%s' is truncated. Continue anyway.\n", fn[heap->i]);

                ks_heapadjust(heap, 0, n, heap);
        }
 
                ks_heapadjust(heap, 0, n, heap);
        }
 

-       for (i = 0; i != n; ++i) {
-               bam_close(fp[i]);
-               free(heap[i].b->data);
-               free(heap[i].b);
-       }
+       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;
        bam_close(fpout);
        free(fp); free(heap);
 }
 int bam_merge(int argc, char *argv[])
 {
        int c, is_by_qname = 0;

-       while ((c = getopt(argc, argv, "n")) >= 0) {
+       char *fn_headers = NULL;
+
+       while ((c = getopt(argc, argv, "h:n")) >= 0) {

                switch (c) {
                switch (c) {

+               case 'h': fn_headers = strdup(optarg); break;

                case 'n': is_by_qname = 1; break;
                }
        }
        if (optind + 2 >= argc) {
                case 'n': is_by_qname = 1; break;
                }
        }
        if (optind + 2 >= argc) {

-               fprintf(stderr, "Usage: samtools merge [-n] <out.bam> <in1.bam> <in2.bam> [...]\n");
+               fprintf(stderr, "\n");
+               fprintf(stderr, "Usage:   samtools merge [-n] [-h inh.sam] <out.bam> <in1.bam> <in2.bam> [...]\n\n");
+               fprintf(stderr, "Options: -n       sort by read 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;
        }
                return 1;
        }

-       bam_merge_core(is_by_qname, argv[optind], argc - optind - 1, argv + optind + 1);
+       bam_merge_core(is_by_qname, argv[optind], fn_headers, argc - optind - 1, argv + optind + 1);
+       free(fn_headers);

        return 0;
 }
 
        return 0;
 }
 


@@ -219,7 +264,7 @@ void bam_sort_core(int is_by_qname, const char *fn, const char *prefix, size_t m
                        fns[i] = (char*)calloc(strlen(prefix) + 20, 1);
                        sprintf(fns[i], "%s.%.4d.bam", prefix, i);
                }
                        fns[i] = (char*)calloc(strlen(prefix) + 20, 1);
                        sprintf(fns[i], "%s.%.4d.bam", prefix, i);
                }

-               bam_merge_core(is_by_qname, fnout, n, fns);
+               bam_merge_core(is_by_qname, fnout, 0, n, fns);

                free(fnout);
                for (i = 0; i < n; ++i) {
                        unlink(fns[i]);
                free(fnout);
                for (i = 0; i < n; ++i) {
                        unlink(fns[i]);

diff --git a/bam_stat.c b/bam_stat.c
index c1c4a432ce8167448f31577d2ce0f650610ff626..ea9deee9d712c71bfd51f7574a5cae055a2c2ace 100644 (file)
--- a/bam_stat.c
+++ b/bam_stat.c
@@ -16,7 +16,7 @@ typedef struct {
                        if ((c)->flag & BAM_FPROPER_PAIR) ++(s)->n_pair_good;           \
                        if ((c)->flag & BAM_FREAD1) ++(s)->n_read1;                                     \
                        if ((c)->flag & BAM_FREAD2) ++(s)->n_read2;                                     \
                        if ((c)->flag & BAM_FPROPER_PAIR) ++(s)->n_pair_good;           \
                        if ((c)->flag & BAM_FREAD1) ++(s)->n_read1;                                     \
                        if ((c)->flag & BAM_FREAD2) ++(s)->n_read2;                                     \

-                       if ((c)->flag & BAM_FMUNMAP) ++(s)->n_sgltn;                            \
+                       if (((c)->flag & BAM_FMUNMAP) && !((c)->flag & BAM_FUNMAP)) ++(s)->n_sgltn;     \

                        if (!((c)->flag & BAM_FUNMAP) && !((c)->flag & BAM_FMUNMAP)) { \
                                ++(s)->n_pair_map;                                                                              \
                                if ((c)->mtid != (c)->tid) {                                                    \
                        if (!((c)->flag & BAM_FUNMAP) && !((c)->flag & BAM_FMUNMAP)) { \
                                ++(s)->n_pair_map;                                                                              \
                                if ((c)->mtid != (c)->tid) {                                                    \

diff --git a/bam_tview.c b/bam_tview.c
index be2579ca0ee28f59792e863bd102b2f7937566d0..4c121e7cbfa012ac0c3fa61d574f03c893e6d271 100644 (file)
--- a/bam_tview.c
+++ b/bam_tview.c
@@ -1,6 +1,21 @@
-#ifndef _NO_CURSES
+#undef _HAVE_CURSES
+
+#if _CURSES_LIB == 0
+#elif _CURSES_LIB == 1

 #include <curses.h>
 #include <curses.h>

-#ifdef NCURSES_VERSION
+#ifndef NCURSES_VERSION
+#warning "_CURSES_LIB=1 but NCURSES_VERSION not defined; tview is NOT compiled"
+#else
+#define _HAVE_CURSES
+#endif
+#elif _CURSES_LIB == 2
+#include <xcurses.h>
+#define _HAVE_CURSES
+#else
+#warning "_CURSES_LIB is not 0, 1 or 2; tview is NOT compiled"
+#endif
+
+#ifdef _HAVE_CURSES

 #include <ctype.h>
 #include <assert.h>
 #include <string.h>
 #include <ctype.h>
 #include <assert.h>
 #include <string.h>


@@ -37,7 +52,7 @@ typedef struct {
        faidx_t *fai;
        bam_maqcns_t *bmc;
 
        faidx_t *fai;
        bam_maqcns_t *bmc;
 

-       int ccol, last_pos, row_shift, base_for, color_for, is_dot, l_ref, ins;
+       int ccol, last_pos, row_shift, base_for, color_for, is_dot, l_ref, ins, no_skip, show_name;

        char *ref;
 } tview_t;
 
        char *ref;
 } tview_t;
 


@@ -50,11 +65,11 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void
        // print referece
        rb = (tv->ref && pos - tv->left_pos < tv->l_ref)? tv->ref[pos - tv->left_pos] : 'N';
        for (i = tv->last_pos + 1; i < pos; ++i) {
        // print referece
        rb = (tv->ref && pos - tv->left_pos < tv->l_ref)? tv->ref[pos - tv->left_pos] : 'N';
        for (i = tv->last_pos + 1; i < pos; ++i) {

-               if (i%10 == 0) mvprintw(0, tv->ccol, "%-d", i+1);
+               if (i%10 == 0 && tv->mcol - tv->ccol >= 10) mvprintw(0, tv->ccol, "%-d", i+1);

                c = tv->ref? tv->ref[i - tv->left_pos] : 'N';
                mvaddch(1, tv->ccol++, c);
        }
                c = tv->ref? tv->ref[i - tv->left_pos] : 'N';
                mvaddch(1, tv->ccol++, c);
        }

-       if (pos%10 == 0) mvprintw(0, tv->ccol, "%-d", pos+1);
+       if (pos%10 == 0 && tv->mcol - tv->ccol >= 10) mvprintw(0, tv->ccol, "%-d", pos+1);

        // print consensus
        call = bam_maqcns_call(n, pl, tv->bmc);
        attr = A_UNDERLINE;
        // print consensus
        call = bam_maqcns_call(n, pl, tv->bmc);
        attr = A_UNDERLINE;


@@ -85,20 +100,28 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void
                                                c = bam_aux_getCSi(p->b, p->qpos);
                                                // assume that if we found one color, we will be able to get the color error
                                                if (tv->is_dot && '-' == bam_aux_getCEi(p->b, p->qpos)) c = bam1_strand(p->b)? ',' : '.';
                                                c = bam_aux_getCSi(p->b, p->qpos);
                                                // assume that if we found one color, we will be able to get the color error
                                                if (tv->is_dot && '-' == bam_aux_getCEi(p->b, p->qpos)) c = bam1_strand(p->b)? ',' : '.';

-                                       }
-                                       else {
-                                               c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)];
-                                               if (tv->is_dot && toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.';
+                                       } else {
+                                               if (tv->show_name) {
+                                                       char *name = bam1_qname(p->b);
+                                                       c = (p->qpos + 1 >= p->b->core.l_qname)? ' ' : name[p->qpos];
+                                               } else {
+                                                       c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)];
+                                                       if (tv->is_dot && toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.';
+                                               }

                                        }
                                } else c = '*';
                        } else { // padding
                                if (j > p->indel) c = '*';
                                else { // insertion
                                        if (tv->base_for ==  TV_BASE_NUCL) {
                                        }
                                } else c = '*';
                        } else { // padding
                                if (j > p->indel) c = '*';
                                else { // insertion
                                        if (tv->base_for ==  TV_BASE_NUCL) {

-                                               c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos + j)];
-                                               if (j == 0 && tv->is_dot && toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.';
-                                       }
-                                       else {
+                                               if (tv->show_name) {
+                                                       char *name = bam1_qname(p->b);
+                                                       c = (p->qpos + j + 1 >= p->b->core.l_qname)? ' ' : name[p->qpos + j];
+                                               } else {
+                                                       c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos + j)];
+                                                       if (j == 0 && tv->is_dot && toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.';
+                                               }
+                                       } else {

                                                c = bam_aux_getCSi(p->b, p->qpos + j);
                                                if (tv->is_dot && '-' == bam_aux_getCEi(p->b, p->qpos + j)) c = bam1_strand(p->b)? ',' : '.';
                                        }
                                                c = bam_aux_getCSi(p->b, p->qpos + j);
                                                if (tv->is_dot && '-' == bam_aux_getCEi(p->b, p->qpos + j)) c = bam1_strand(p->b)? ',' : '.';
                                        }


@@ -177,13 +200,10 @@ tview_t *tv_init(const char *fn, const char *fn_fa)
        clear();
        noecho();
        cbreak();
        clear();
        noecho();
        cbreak();

-#ifdef NCURSES_VERSION
+       tv->mrow = 24; tv->mcol = 80;

        getmaxyx(stdscr, tv->mrow, tv->mcol);
        getmaxyx(stdscr, tv->mrow, tv->mcol);

-#else
-       tv->mrow = 80; tv->mcol = 40;
-#endif

        tv->wgoto = newwin(3, TV_MAX_GOTO + 10, 10, 5);
        tv->wgoto = newwin(3, TV_MAX_GOTO + 10, 10, 5);

-       tv->whelp = newwin(27, 40, 5, 5);
+       tv->whelp = newwin(29, 40, 5, 5);

        tv->color_for = TV_COLOR_MAPQ;
        start_color();
        init_pair(1, COLOR_BLUE, COLOR_BLACK);
        tv->color_for = TV_COLOR_MAPQ;
        start_color();
        init_pair(1, COLOR_BLUE, COLOR_BLACK);


@@ -216,6 +236,14 @@ void tv_destroy(tview_t *tv)
 int tv_fetch_func(const bam1_t *b, void *data)
 {
        tview_t *tv = (tview_t*)data;
 int tv_fetch_func(const bam1_t *b, void *data)
 {
        tview_t *tv = (tview_t*)data;

+       if (tv->no_skip) {
+               uint32_t *cigar = bam1_cigar(b); // this is cheating...
+               int i;
+               for (i = 0; i <b->core.n_cigar; ++i) {
+                       if ((cigar[i]&0xf) == BAM_CREF_SKIP)
+                               cigar[i] = cigar[i]>>4<<4 | BAM_CDEL;
+               }
+       }

        bam_lplbuf_push(b, tv->lplbuf);
        return 0;
 }
        bam_lplbuf_push(b, tv->lplbuf);
        return 0;
 }


@@ -240,6 +268,13 @@ int tv_draw_aln(tview_t *tv, int tid, int pos)
        bam_lplbuf_reset(tv->lplbuf);
        bam_fetch(tv->fp, tv->idx, tv->curr_tid, tv->left_pos, tv->left_pos + tv->mcol, tv, tv_fetch_func);
        bam_lplbuf_push(0, tv->lplbuf);
        bam_lplbuf_reset(tv->lplbuf);
        bam_fetch(tv->fp, tv->idx, tv->curr_tid, tv->left_pos, tv->left_pos + tv->mcol, tv, tv_fetch_func);
        bam_lplbuf_push(0, tv->lplbuf);

+
+       while (tv->ccol < tv->mcol) {
+               int pos = tv->last_pos + 1;
+               if (pos%10 == 0 && tv->mcol - tv->ccol >= 10) mvprintw(0, tv->ccol, "%-d", pos+1);
+               mvaddch(1, tv->ccol++, (tv->ref && pos < tv->l_ref)? tv->ref[pos - tv->left_pos] : 'N');
+               ++tv->last_pos;
+       }

        return 0;
 }
 
        return 0;
 }
 


@@ -292,6 +327,8 @@ static void tv_win_help(tview_t *tv) {
        mvwprintw(win, r++, 2, "c          Color for cs color");
        mvwprintw(win, r++, 2, "z          Color for cs qual");
        mvwprintw(win, r++, 2, ".          Toggle on/off dot view");
        mvwprintw(win, r++, 2, "c          Color for cs color");
        mvwprintw(win, r++, 2, "z          Color for cs qual");
        mvwprintw(win, r++, 2, ".          Toggle on/off dot view");

+       mvwprintw(win, r++, 2, "s          Toggle on/off ref skip");
+       mvwprintw(win, r++, 2, "r          Toggle on/off rd name");

        mvwprintw(win, r++, 2, "N          Turn on nt view");
        mvwprintw(win, r++, 2, "C          Turn on cs view");
        mvwprintw(win, r++, 2, "i          Toggle on/off ins");
        mvwprintw(win, r++, 2, "N          Turn on nt view");
        mvwprintw(win, r++, 2, "C          Turn on cs view");
        mvwprintw(win, r++, 2, "i          Toggle on/off ins");


@@ -310,7 +347,6 @@ void tv_loop(tview_t *tv)
        tid = tv->curr_tid; pos = tv->left_pos;
        while (1) {
                int c = getch();
        tid = tv->curr_tid; pos = tv->left_pos;
        while (1) {
                int c = getch();

-               //if(256 < c) {c = 1 + (c%256);} // Terminal was displaying ctrl-H as 263 via ssh from Mac OS X 10.5 computer 

                switch (c) {
                        case '?': tv_win_help(tv); break;
                        case '\033':
                switch (c) {
                        case '?': tv_win_help(tv); break;
                        case '\033':


@@ -321,6 +357,8 @@ void tv_loop(tview_t *tv)
                        case 'n': tv->color_for = TV_COLOR_NUCL; break;
                        case 'c': tv->color_for = TV_COLOR_COL; break;
                        case 'z': tv->color_for = TV_COLOR_COLQ; break;
                        case 'n': tv->color_for = TV_COLOR_NUCL; break;
                        case 'c': tv->color_for = TV_COLOR_COL; break;
                        case 'z': tv->color_for = TV_COLOR_COLQ; break;

+                       case 's': tv->no_skip = !tv->no_skip; break;
+                       case 'r': tv->show_name = !tv->show_name; break;

                        case KEY_LEFT:
                        case 'h': --pos; break;
                        case KEY_RIGHT:
                        case KEY_LEFT:
                        case 'h': --pos; break;
                        case KEY_RIGHT:


@@ -342,9 +380,7 @@ void tv_loop(tview_t *tv)
                        case 'k': ++tv->row_shift; break;
                        case KEY_BACKSPACE:
                        case '\177': pos -= tv->mcol; break;
                        case 'k': ++tv->row_shift; break;
                        case KEY_BACKSPACE:
                        case '\177': pos -= tv->mcol; break;

-#ifdef KEY_RESIZE

                        case KEY_RESIZE: getmaxyx(stdscr, tv->mrow, tv->mcol); break;
                        case KEY_RESIZE: getmaxyx(stdscr, tv->mrow, tv->mcol); break;

-#endif

                        default: continue;
                }
                if (pos < 0) pos = 0;
                        default: continue;
                }
                if (pos < 0) pos = 0;


@@ -368,12 +404,12 @@ int bam_tview_main(int argc, char *argv[])
        tv_destroy(tv);
        return 0;
 }
        tv_destroy(tv);
        return 0;
 }

-#else // #ifdef NCURSES_VERSION
-#warning "The ncurses library is unavailable; tview is disabled."
+#else // #ifdef _HAVE_CURSES
+#include <stdio.h>
+#warning "No curses library is available; tview is disabled."

 int bam_tview_main(int argc, char *argv[])
 {
        fprintf(stderr, "[bam_tview_main] The ncurses library is unavailable; tview is not compiled.\n");
        return 1;
 }
 int bam_tview_main(int argc, char *argv[])
 {
        fprintf(stderr, "[bam_tview_main] The ncurses library is unavailable; tview is not compiled.\n");
        return 1;
 }

-#endif
-#endif // #ifndef _NO_CURSES
+#endif // #ifdef _HAVE_CURSES

diff --git a/bamtk.c b/bamtk.c
index 3386836e04e3a0b44eb47ff9d925d0db20b8ea90..ea6667205c9b642e9e5540470d33b5ab274d727a 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -1,10 +1,15 @@
 #include <stdio.h>
 #include <unistd.h>
 #include <assert.h>
 #include <stdio.h>
 #include <unistd.h>
 #include <assert.h>

+#include <fcntl.h>

 #include "bam.h"
 
 #include "bam.h"
 

+#ifdef _USE_KNETFILE
+#include "knetfile.h"
+#endif
+

 #ifndef PACKAGE_VERSION
 #ifndef PACKAGE_VERSION

-#define PACKAGE_VERSION "0.1.5c (r385)"
+#define PACKAGE_VERSION "0.1.6 (r453)"

 #endif
 
 int bam_taf2baf(int argc, char *argv[]);
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);


@@ -71,27 +76,33 @@ static int usage()
        fprintf(stderr, "Program: samtools (Tools for alignments in the SAM format)\n");
        fprintf(stderr, "Version: %s\n\n", PACKAGE_VERSION);
        fprintf(stderr, "Usage:   samtools <command> [options]\n\n");
        fprintf(stderr, "Program: samtools (Tools for alignments in the SAM format)\n");
        fprintf(stderr, "Version: %s\n\n", PACKAGE_VERSION);
        fprintf(stderr, "Usage:   samtools <command> [options]\n\n");

-       fprintf(stderr, "Command: import      import from SAM (obsolete; use `view')\n");
-       fprintf(stderr, "         view        export to the text format\n");
+       fprintf(stderr, "Command: view        SAM<->BAM conversion\n");

        fprintf(stderr, "         sort        sort alignment file\n");
        fprintf(stderr, "         sort        sort alignment file\n");

-       fprintf(stderr, "         merge       merge multiple sorted alignment files\n");

        fprintf(stderr, "         pileup      generate pileup output\n");
        fprintf(stderr, "         faidx       index/extract FASTA\n");
        fprintf(stderr, "         pileup      generate pileup output\n");
        fprintf(stderr, "         faidx       index/extract FASTA\n");

-#ifndef _NO_CURSES
+#if _CURSES_LIB != 0

        fprintf(stderr, "         tview       text alignment viewer\n");
 #endif
        fprintf(stderr, "         index       index alignment\n");
        fprintf(stderr, "         fixmate     fix mate information\n");
        fprintf(stderr, "         tview       text alignment viewer\n");
 #endif
        fprintf(stderr, "         index       index alignment\n");
        fprintf(stderr, "         fixmate     fix mate information\n");

-       fprintf(stderr, "         rmdup       remove PCR duplicates\n");

        fprintf(stderr, "         glfview     print GLFv3 file\n");
        fprintf(stderr, "         flagstat    simple stats\n");
        fprintf(stderr, "         glfview     print GLFv3 file\n");
        fprintf(stderr, "         flagstat    simple stats\n");

-       fprintf(stderr, "         fillmd      fill the MD tag and change identical base to =\n");
+       fprintf(stderr, "         calmd       recalculate MD/NM tags and '=' bases\n");
+       fprintf(stderr, "         merge       merge sorted alignments (Picard recommended)\n");
+       fprintf(stderr, "         rmdup       remove PCR duplicates (Picard recommended)\n");

        fprintf(stderr, "\n");
        return 1;
 }
 
 int main(int argc, char *argv[])
 {
        fprintf(stderr, "\n");
        return 1;
 }
 
 int main(int argc, char *argv[])
 {

+#ifdef _WIN32
+       setmode(fileno(stdout), O_BINARY);
+       setmode(fileno(stdin),  O_BINARY);
+#ifdef _USE_KNETFILE
+       knet_win32_init();
+#endif
+#endif

        if (argc < 2) return usage();
        if (strcmp(argv[1], "view") == 0) return main_samview(argc-1, argv+1);
        else if (strcmp(argv[1], "import") == 0) return main_import(argc-1, argv+1);
        if (argc < 2) return usage();
        if (strcmp(argv[1], "view") == 0) return main_samview(argc-1, argv+1);
        else if (strcmp(argv[1], "import") == 0) return main_import(argc-1, argv+1);


@@ -106,8 +117,9 @@ int main(int argc, char *argv[])
        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], "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);
        else if (strcmp(argv[1], "fillmd") == 0) return bam_fillmd(argc-1, argv+1);

-#ifndef _NO_CURSES
+#if _CURSES_LIB != 0

        else if (strcmp(argv[1], "tview") == 0) return bam_tview_main(argc-1, argv+1);
 #endif
        else {
        else if (strcmp(argv[1], "tview") == 0) return bam_tview_main(argc-1, argv+1);
 #endif
        else {

diff --git a/bgzf.c b/bgzf.c
index fe4e31da109f6ca523e73088fb3fd27aa9672587..646b2b4a142e2c35f4ab039b515e68fcd9126d66 100644 (file)
--- a/bgzf.c
+++ b/bgzf.c
@@ -1,13 +1,25 @@
-/*
- * The Broad Institute
- * SOFTWARE COPYRIGHT NOTICE AGREEMENT
- * This software and its documentation are copyright 2008 by the
- * Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
- *
- * This software is supplied without any warranty or guaranteed support whatsoever.
- * Neither the Broad Institute nor MIT can be responsible for its use, misuse,
- * or functionality.
- */
+/* The MIT License
+
+   Copyright (c) 2008 Broad Institute / Massachusetts Institute of Technology
+
+   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.
+*/

 
 /*
   2009-06-29 by lh3: cache recent uncompressed blocks.
 
 /*
   2009-06-29 by lh3: cache recent uncompressed blocks.


@@ -31,10 +43,15 @@ typedef struct {
 } cache_t;
 KHASH_MAP_INIT_INT64(cache, cache_t)
 
 } cache_t;
 KHASH_MAP_INIT_INT64(cache, cache_t)
 

+#if defined(_WIN32) || defined(_MSC_VER)
+#define ftello(fp) ftell(fp)
+#define fseeko(fp, offset, whence) fseek(fp, offset, whence)
+#else

 extern off_t ftello(FILE *stream);
 extern int fseeko(FILE *stream, off_t offset, int whence);
 extern off_t ftello(FILE *stream);
 extern int fseeko(FILE *stream, off_t offset, int whence);

+#endif

 
 

-typedef int8_t byte;
+typedef int8_t bgzf_byte_t;

 
 static const int DEFAULT_BLOCK_SIZE = 64 * 1024;
 static const int MAX_BLOCK_SIZE = 64 * 1024;
 
 static const int DEFAULT_BLOCK_SIZE = 64 * 1024;
 static const int MAX_BLOCK_SIZE = 64 * 1024;


@@ -81,9 +98,9 @@ packInt32(uint8_t* buffer, uint32_t value)
     buffer[3] = value >> 24;
 }
 
     buffer[3] = value >> 24;
 }
 

-inline
+static inline

 int
 int

-min(int x, int y)
+bgzf_min(int x, int y)

 {
     return (x < y) ? x : y;
 }
 {
     return (x < y) ? x : y;
 }


@@ -169,14 +186,20 @@ bgzf_open(const char* __restrict path, const char* __restrict mode)
                fp->open_mode = 'r';
                fp->x.fpr = file;
 #else
                fp->open_mode = 'r';
                fp->x.fpr = file;
 #else

-               int oflag = O_RDONLY;
-               int fd = open(path, oflag);
+               int fd, oflag = O_RDONLY;
+#ifdef _WIN32
+               oflag |= O_BINARY;
+#endif
+               fd = open(path, oflag);

                if (fd == -1) return 0;
         fp = open_read(fd);
 #endif
     } else if (mode[0] == 'w' || mode[0] == 'W') {
                if (fd == -1) return 0;
         fp = open_read(fd);
 #endif
     } else if (mode[0] == 'w' || mode[0] == 'W') {

-               int oflag = O_WRONLY | O_CREAT | O_TRUNC;
-               int fd = open(path, oflag, 0644);
+               int fd, oflag = O_WRONLY | O_CREAT | O_TRUNC;
+#ifdef _WIN32
+               oflag |= O_BINARY;
+#endif
+               fd = open(path, oflag, 0644);

                if (fd == -1) return 0;
         fp = open_write(fd, strstr(mode, "u")? 1 : 0);
     }
                if (fd == -1) return 0;
         fp = open_write(fd, strstr(mode, "u")? 1 : 0);
     }


@@ -206,7 +229,7 @@ deflate_block(BGZF* fp, int block_length)
     // Deflate the block in fp->uncompressed_block into fp->compressed_block.
     // Also adds an extra field that stores the compressed block length.
 
     // Deflate the block in fp->uncompressed_block into fp->compressed_block.
     // Also adds an extra field that stores the compressed block length.
 

-    byte* buffer = fp->compressed_block;
+    bgzf_byte_t* buffer = fp->compressed_block;

     int buffer_size = fp->compressed_block_size;
 
     // Init gzip header
     int buffer_size = fp->compressed_block_size;
 
     // Init gzip header


@@ -337,10 +360,10 @@ inflate_block(BGZF* fp, int block_length)
 
 static
 int
 
 static
 int

-check_header(const byte* header)
+check_header(const bgzf_byte_t* header)

 {
     return (header[0] == GZIP_ID1 &&
 {
     return (header[0] == GZIP_ID1 &&

-            header[1] == (byte) GZIP_ID2 &&
+            header[1] == (bgzf_byte_t) GZIP_ID2 &&

             header[2] == Z_DEFLATED &&
             (header[3] & FLG_FEXTRA) != 0 &&
             unpackInt16((uint8_t*)&header[10]) == BGZF_XLEN &&
             header[2] == Z_DEFLATED &&
             (header[3] & FLG_FEXTRA) != 0 &&
             unpackInt16((uint8_t*)&header[10]) == BGZF_XLEN &&


@@ -410,7 +433,7 @@ static
 int
 read_block(BGZF* fp)
 {
 int
 read_block(BGZF* fp)
 {

-    byte header[BLOCK_HEADER_LENGTH];
+    bgzf_byte_t header[BLOCK_HEADER_LENGTH];

        int size = 0;
 #ifdef _USE_KNETFILE
     int64_t block_address = knet_tell(fp->x.fpr);
        int size = 0;
 #ifdef _USE_KNETFILE
     int64_t block_address = knet_tell(fp->x.fpr);


@@ -435,7 +458,7 @@ read_block(BGZF* fp)
         return -1;
     }
     int block_length = unpackInt16((uint8_t*)&header[16]) + 1;
         return -1;
     }
     int block_length = unpackInt16((uint8_t*)&header[16]) + 1;

-    byte* compressed_block = (byte*) fp->compressed_block;
+    bgzf_byte_t* compressed_block = (bgzf_byte_t*) fp->compressed_block;

     memcpy(compressed_block, header, BLOCK_HEADER_LENGTH);
     int remaining = block_length - BLOCK_HEADER_LENGTH;
 #ifdef _USE_KNETFILE
     memcpy(compressed_block, header, BLOCK_HEADER_LENGTH);
     int remaining = block_length - BLOCK_HEADER_LENGTH;
 #ifdef _USE_KNETFILE


@@ -474,7 +497,7 @@ bgzf_read(BGZF* fp, void* data, int length)
     }
 
     int bytes_read = 0;
     }
 
     int bytes_read = 0;

-    byte* output = data;
+    bgzf_byte_t* output = data;

     while (bytes_read < length) {
         int available = fp->block_length - fp->block_offset;
         if (available <= 0) {
     while (bytes_read < length) {
         int available = fp->block_length - fp->block_offset;
         if (available <= 0) {


@@ -486,8 +509,8 @@ bgzf_read(BGZF* fp, void* data, int length)
                 break;
             }
         }
                 break;
             }
         }

-        int copy_length = min(length-bytes_read, available);
-        byte* buffer = fp->uncompressed_block;
+        int copy_length = bgzf_min(length-bytes_read, available);
+        bgzf_byte_t* buffer = fp->uncompressed_block;

         memcpy(output, buffer + fp->block_offset, copy_length);
         fp->block_offset += copy_length;
         output += copy_length;
         memcpy(output, buffer + fp->block_offset, copy_length);
         fp->block_offset += copy_length;
         output += copy_length;


@@ -540,12 +563,12 @@ bgzf_write(BGZF* fp, const void* data, int length)
         fp->uncompressed_block = malloc(fp->uncompressed_block_size);
     }
 
         fp->uncompressed_block = malloc(fp->uncompressed_block_size);
     }
 

-    const byte* input = data;
+    const bgzf_byte_t* input = data;

     int block_length = fp->uncompressed_block_size;
     int bytes_written = 0;
     while (bytes_written < length) {
     int block_length = fp->uncompressed_block_size;
     int bytes_written = 0;
     while (bytes_written < length) {

-        int copy_length = min(block_length - fp->block_offset, length - bytes_written);
-        byte* buffer = fp->uncompressed_block;
+        int copy_length = bgzf_min(block_length - fp->block_offset, length - bytes_written);
+        bgzf_byte_t* buffer = fp->uncompressed_block;

         memcpy(buffer + fp->block_offset, input, copy_length);
         fp->block_offset += copy_length;
         input += copy_length;
         memcpy(buffer + fp->block_offset, input, copy_length);
         fp->block_offset += copy_length;
         input += copy_length;


@@ -566,6 +589,14 @@ bgzf_close(BGZF* fp)
         if (flush_block(fp) != 0) {
             return -1;
         }
         if (flush_block(fp) != 0) {
             return -1;
         }

+               { // add an empty block
+                       int count, block_length = deflate_block(fp, 0);
+#ifdef _USE_KNETFILE
+                       count = fwrite(fp->compressed_block, 1, block_length, fp->x.fpw);
+#else
+                       count = fwrite(fp->compressed_block, 1, block_length, fp->file);
+#endif
+               }

 #ifdef _USE_KNETFILE
         if (fflush(fp->x.fpw) != 0) {
 #else
 #ifdef _USE_KNETFILE
         if (fflush(fp->x.fpw) != 0) {
 #else


@@ -605,6 +636,25 @@ void bgzf_set_cache_size(BGZF *fp, int cache_size)
        if (fp) fp->cache_size = cache_size;
 }
 
        if (fp) fp->cache_size = cache_size;
 }
 

+int bgzf_check_EOF(BGZF *fp)
+{
+       static uint8_t magic[28] = "\037\213\010\4\0\0\0\0\0\377\6\0\102\103\2\0\033\0\3\0\0\0\0\0\0\0\0\0";
+       uint8_t buf[28];
+       off_t offset;
+#ifdef _USE_KNETFILE
+       offset = knet_tell(fp->x.fpr);
+       if (knet_seek(fp->x.fpr, -28, SEEK_END) != 0) return -1;
+       knet_read(fp->x.fpr, buf, 28);
+       knet_seek(fp->x.fpr, offset, SEEK_SET);
+#else
+       offset = ftello(fp->file);
+       if (fseeko(fp->file, -28, SEEK_END) != 0) return -1;
+       fread(buf, 1, 28, fp->file);
+       fseeko(fp->file, offset, SEEK_SET);
+#endif
+       return (memcmp(magic, buf, 28) == 0)? 1 : 0;
+}
+

 int64_t
 bgzf_seek(BGZF* fp, int64_t pos, int where)
 {
 int64_t
 bgzf_seek(BGZF* fp, int64_t pos, int where)
 {


@@ -631,4 +681,3 @@ bgzf_seek(BGZF* fp, int64_t pos, int where)
     fp->block_offset = block_offset;
     return 0;
 }
     fp->block_offset = block_offset;
     return 0;
 }

-

diff --git a/bgzf.h b/bgzf.h
index d5eeafe7e527c57823b8b68cb2c32511f1b0fadd..91b33177ff07c667b16f4db2ab724081beee5499 100644 (file)
--- a/bgzf.h
+++ b/bgzf.h
@@ -1,13 +1,25 @@
-/*
- * The Broad Institute
- * SOFTWARE COPYRIGHT NOTICE AGREEMENT
- * This software and its documentation are copyright 2008 by the
- * Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
- *
- * This software is supplied without any warranty or guaranteed support whatsoever.
- * Neither the Broad Institute nor MIT can be responsible for its use, misuse,
- * or functionality.
- */
+/* The MIT License
+
+   Copyright (c) 2008 Broad Institute / Massachusetts Institute of Technology
+
+   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 __BGZF_H
 #define __BGZF_H
 
 #ifndef __BGZF_H
 #define __BGZF_H


@@ -113,6 +125,8 @@ int64_t bgzf_seek(BGZF* fp, int64_t pos, int where);
  */
 void bgzf_set_cache_size(BGZF *fp, int cache_size);
 
  */
 void bgzf_set_cache_size(BGZF *fp, int cache_size);
 

+int bgzf_check_EOF(BGZF *fp);
+

 #ifdef __cplusplus
 }
 #endif
 #ifdef __cplusplus
 }
 #endif

diff --git a/bgzip.c b/bgzip.c
index c58d55d193c7daf7378f6a763ee5664f6f71461e..eb88195ec279a99f63b9a6780103d00f7e5829ca 100644 (file)
--- a/bgzip.c
+++ b/bgzip.c
@@ -1,13 +1,26 @@
-/*
- * The Broad Institute
- * SOFTWARE COPYRIGHT NOTICE AGREEMENT
- * This software and its documentation are copyright 2008 by the
- * Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
- *
- * This software is supplied without any warranty or guaranteed support whatsoever.
- * Neither the Broad Institute nor MIT can be responsible for its use, misuse,
- * or functionality.
- */
+/* The MIT License
+
+   Copyright (c) 2008 Broad Institute / Massachusetts Institute of Technology
+
+   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 <string.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <stdio.h>

diff --git a/faidx.c b/faidx.c
index 36366c20db595ee9ffeaac4d87628dc40c3690d2..055445f266cbdb429fe65c2aeb58c763af86d1d9 100644 (file)
--- a/faidx.c
+++ b/faidx.c
@@ -14,8 +14,13 @@ KHASH_MAP_INIT_STR(s, faidx1_t)
 #ifndef _NO_RAZF
 #include "razf.h"
 #else
 #ifndef _NO_RAZF
 #include "razf.h"
 #else

+#ifdef _WIN32
+#define ftello(fp) ftell(fp)
+#define fseeko(fp, offset, whence) fseek(fp, offset, whence)
+#else

 extern off_t ftello(FILE *stream);
 extern int fseeko(FILE *stream, off_t offset, int whence);
 extern off_t ftello(FILE *stream);
 extern int fseeko(FILE *stream, off_t offset, int whence);

+#endif

 #define RAZF FILE
 #define razf_read(fp, buf, size) fread(buf, 1, size, fp)
 #define razf_open(fn, mode) fopen(fn, mode)
 #define RAZF FILE
 #define razf_read(fp, buf, size) fread(buf, 1, size, fp)
 #define razf_open(fn, mode) fopen(fn, mode)


@@ -134,7 +139,11 @@ void fai_save(const faidx_t *fai, FILE *fp)
                faidx1_t x;
                k = kh_get(s, fai->hash, fai->name[i]);
                x = kh_value(fai->hash, k);
                faidx1_t x;
                k = kh_get(s, fai->hash, fai->name[i]);
                x = kh_value(fai->hash, k);

+#ifdef _WIN32
+               fprintf(fp, "%s\t%d\t%ld\t%d\t%d\n", fai->name[i], (int)x.len, (long)x.offset, (int)x.line_blen, (int)x.line_len);
+#else

                fprintf(fp, "%s\t%d\t%lld\t%d\t%d\n", fai->name[i], (int)x.len, (long long)x.offset, (int)x.line_blen, (int)x.line_len);
                fprintf(fp, "%s\t%d\t%lld\t%d\t%d\n", fai->name[i], (int)x.len, (long long)x.offset, (int)x.line_blen, (int)x.line_len);

+#endif

        }
 }
 
        }
 }
 


@@ -143,14 +152,22 @@ faidx_t *fai_read(FILE *fp)
        faidx_t *fai;
        char *buf, *p;
        int len, line_len, line_blen;
        faidx_t *fai;
        char *buf, *p;
        int len, line_len, line_blen;

+#ifdef _WIN32
+       long offset;
+#else

        long long offset;
        long long offset;

+#endif

        fai = (faidx_t*)calloc(1, sizeof(faidx_t));
        fai->hash = kh_init(s);
        buf = (char*)calloc(0x10000, 1);
        while (!feof(fp) && fgets(buf, 0x10000, fp)) {
                for (p = buf; *p && isgraph(*p); ++p);
                *p = 0; ++p;
        fai = (faidx_t*)calloc(1, sizeof(faidx_t));
        fai->hash = kh_init(s);
        buf = (char*)calloc(0x10000, 1);
        while (!feof(fp) && fgets(buf, 0x10000, fp)) {
                for (p = buf; *p && isgraph(*p); ++p);
                *p = 0; ++p;

+#ifdef _WIN32
+               sscanf(p, "%d%ld%d%d", &len, &offset, &line_blen, &line_len);
+#else

                sscanf(p, "%d%lld%d%d", &len, &offset, &line_blen, &line_len);
                sscanf(p, "%d%lld%d%d", &len, &offset, &line_blen, &line_len);

+#endif

                fai_insert_index(fai, buf, len, line_len, line_blen, offset);
        }
        free(buf);
                fai_insert_index(fai, buf, len, line_len, line_blen, offset);
        }
        free(buf);


@@ -183,7 +200,7 @@ int fai_build(const char *fn)
        }
        fai = fai_build_core(rz);
        razf_close(rz);
        }
        fai = fai_build_core(rz);
        razf_close(rz);

-       fp = fopen(str, "w");
+       fp = fopen(str, "wb");

        if (fp == 0) {
                fprintf(stderr, "[fai_build] fail to write FASTA index.\n");
                fai_destroy(fai); free(str);
        if (fp == 0) {
                fprintf(stderr, "[fai_build] fail to write FASTA index.\n");
                fai_destroy(fai); free(str);


@@ -203,7 +220,7 @@ faidx_t *fai_load(const char *fn)
        faidx_t *fai;
        str = (char*)calloc(strlen(fn) + 5, 1);
        sprintf(str, "%s.fai", fn);
        faidx_t *fai;
        str = (char*)calloc(strlen(fn) + 5, 1);
        sprintf(str, "%s.fai", fn);

-       fp = fopen(str, "r");
+       fp = fopen(str, "rb");

        if (fp == 0) {
                fprintf(stderr, "[fai_load] build FASTA index.\n");
                fai_build(fn);
        if (fp == 0) {
                fprintf(stderr, "[fai_load] build FASTA index.\n");
                fai_build(fn);


@@ -216,7 +233,7 @@ faidx_t *fai_load(const char *fn)
        }
        fai = fai_read(fp);
        fclose(fp);
        }
        fai = fai_read(fp);
        fclose(fp);

-       fai->rz = razf_open(fn, "r");
+       fai->rz = razf_open(fn, "rb");

        free(str);
        if (fai->rz == 0) {
                fprintf(stderr, "[fai_load] fail to open FASTA file.\n");
        free(str);
        if (fai->rz == 0) {
                fprintf(stderr, "[fai_load] fail to open FASTA file.\n");

diff --git a/knetfile.c b/knetfile.c
index cef197dd08b0b29e97be7a3dade6cac54d244c46..e110aa720f6577e7f8dee34fddc32cb6220d7154 100644 (file)
--- a/knetfile.c
+++ b/knetfile.c
@@ -1,15 +1,64 @@
+/* The MIT License
+
+   Copyright (c) 2008 Genome Research Ltd (GRL).
+
+   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.
+*/
+
+/* Contact: Heng Li <lh3@sanger.ac.uk> */
+
+/* Probably I will not do socket programming in the next few years and
+   therefore I decide to heavily annotate this file, for Linux and
+   Windows as well.  -lh3 */
+

 #include <time.h>
 #include <stdio.h>
 #include <time.h>
 #include <stdio.h>

-#include <netdb.h>

 #include <ctype.h>
 #include <stdlib.h>
 #include <string.h>
 #include <unistd.h>
 #include <sys/types.h>
 #include <ctype.h>
 #include <stdlib.h>
 #include <string.h>
 #include <unistd.h>
 #include <sys/types.h>

+
+#ifdef _WIN32
+#include <winsock.h>
+#else
+#include <netdb.h>

 #include <arpa/inet.h>
 #include <sys/socket.h>
 #include <arpa/inet.h>
 #include <sys/socket.h>

+#endif
+

 #include "knetfile.h"
 
 #include "knetfile.h"
 

+/* In winsock.h, the type of a socket is SOCKET, which is: "typedef
+ * u_int SOCKET". An invalid SOCKET is: "(SOCKET)(~0)", or signed
+ * integer -1. In knetfile.c, I use "int" for socket type
+ * throughout. This should be improved to avoid confusion.
+ *
+ * In Linux/Mac, recv() and read() do almost the same thing. You can see
+ * in the header file that netread() is simply an alias of read(). In
+ * Windows, however, they are different and using recv() is mandatory.
+ */
+
+/* This function tests if the file handler is ready for reading (or
+ * writing if is_read==0). */

 static int socket_wait(int fd, int is_read)
 {
        fd_set fds, *fdr = 0, *fdw = 0;
 static int socket_wait(int fd, int is_read)
 {
        fd_set fds, *fdr = 0, *fdw = 0;


@@ -25,13 +74,119 @@ static int socket_wait(int fd, int is_read)
        return ret;
 }
 
        return ret;
 }
 

+#ifndef _WIN32
+/* This function does not work with Windows due to the lack of
+ * getaddrinfo() in winsock. It is addapted from an example in "Beej's
+ * Guide to Network Programming" (http://beej.us/guide/bgnet/). */
+static int socket_connect(const char *host, const char *port)
+{
+#define __err_connect(func) do { perror(func); freeaddrinfo(res); return -1; } while (0)
+
+       int on = 1, fd;
+       struct linger lng = { 0, 0 };
+       struct addrinfo hints, *res;
+       memset(&hints, 0, sizeof(struct addrinfo));
+       hints.ai_family = AF_UNSPEC;
+       hints.ai_socktype = SOCK_STREAM;
+       /* In Unix/Mac, getaddrinfo() is the most convenient way to get
+        * server information. */
+       if (getaddrinfo(host, port, &hints, &res) != 0) __err_connect("getaddrinfo");
+       if ((fd = socket(res->ai_family, res->ai_socktype, res->ai_protocol)) == -1) __err_connect("socket");
+       /* The following two setsockopt() are used by ftplib
+        * (http://nbpfaus.net/~pfau/ftplib/). I am not sure if they
+        * necessary. */
+       if (setsockopt(fd, SOL_SOCKET, SO_REUSEADDR, &on, sizeof(on)) == -1) __err_connect("setsockopt");
+       if (setsockopt(fd, SOL_SOCKET, SO_LINGER, &lng, sizeof(lng)) == -1) __err_connect("setsockopt");
+       if (connect(fd, res->ai_addr, res->ai_addrlen) != 0) __err_connect("connect");
+       freeaddrinfo(res);
+       return fd;
+}
+#else
+/* MinGW's printf has problem with "%lld" */
+char *uint64tostr(char *buf, uint64_t x)
+{
+       int i, cnt;
+       for (i = 0; x; x /= 10) buf[i++] = '0' + x%10;
+       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;
+}
+/* In windows, the first thing is to establish the TCP connection. */
+int knet_win32_init()
+{
+       WSADATA wsaData;
+       return WSAStartup(MAKEWORD(2, 2), &wsaData);
+}
+void knet_win32_destroy()
+{
+       WSACleanup();
+}
+/* A slightly modfied version of the following function also works on
+ * Mac (and presummably Linux). However, this function is not stable on
+ * my Mac. It sometimes works fine but sometimes does not. Therefore for
+ * 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)
+
+       int on = 1;
+       SOCKET fd;
+       struct linger lng = { 0, 0 };
+       struct sockaddr_in server;
+       struct hostent *hp = 0;
+       // open socket
+       if ((fd = socket(AF_INET, SOCK_STREAM, IPPROTO_TCP)) == INVALID_SOCKET) __err_connect("socket");
+       if (setsockopt(fd, SOL_SOCKET, SO_REUSEADDR, (char*)&on, sizeof(on)) == -1) __err_connect("setsockopt");
+       if (setsockopt(fd, SOL_SOCKET, SO_LINGER, (char*)&lng, sizeof(lng)) == -1) __err_connect("setsockopt");
+       // get host info
+       if (isalpha(host[0])) hp = gethostbyname(host);
+       else {
+               struct in_addr addr;
+               addr.s_addr = inet_addr(host);
+               hp = gethostbyaddr((char*)&addr, 4, AF_INET);
+       }
+       if (hp == 0) __err_connect("gethost");
+       // connect
+       server.sin_addr.s_addr = *((unsigned long*)hp->h_addr);
+       server.sin_family= AF_INET;
+       server.sin_port = htons(atoi(port));
+       if (connect(fd, (struct sockaddr*)&server, sizeof(server)) != 0) __err_connect("connect");
+       // freehostent(hp); // strangely in MSDN, hp is NOT freed (memory leak?!)
+       return fd;
+}
+#endif
+
+static off_t my_netread(int fd, void *buf, off_t len)
+{
+       off_t rest = len, curr, l = 0;
+       /* recv() and read() may not read the required length of data with
+        * one call. They have to be called repeatedly. */
+       while (rest) {
+               if (socket_wait(fd, 1) <= 0) break; // socket is not ready for reading
+               curr = netread(fd, buf + l, rest);
+               /* According to the glibc manual, section 13.2, a zero returned
+                * value indicates end-of-file (EOF), which should mean that
+                * read() will not return zero if EOF has not been met but data
+                * are not immediately available. */
+               if (curr == 0) break;
+               l += curr; rest -= curr;
+       }
+       return l;
+}
+
+/*************************
+ * FTP specific routines *
+ *************************/
+

 static int kftp_get_response(knetFile *ftp)
 {
        unsigned char c;
        int n = 0;
        char *p;
        if (socket_wait(ftp->ctrl_fd, 1) <= 0) return 0;
 static int kftp_get_response(knetFile *ftp)
 {
        unsigned char c;
        int n = 0;
        char *p;
        if (socket_wait(ftp->ctrl_fd, 1) <= 0) return 0;

-       while (read(ftp->ctrl_fd, &c, 1)) { // FIXME: this is *VERY BAD* for unbuffered I/O
+       while (netread(ftp->ctrl_fd, &c, 1)) { // FIXME: this is *VERY BAD* for unbuffered I/O

                //fputc(c, stderr);
                if (n >= ftp->max_response) {
                        ftp->max_response = ftp->max_response? ftp->max_response<<1 : 256;
                //fputc(c, stderr);
                if (n >= ftp->max_response) {
                        ftp->max_response = ftp->max_response? ftp->max_response<<1 : 256;


@@ -53,7 +208,7 @@ static int kftp_get_response(knetFile *ftp)
 static int kftp_send_cmd(knetFile *ftp, const char *cmd, int is_get)
 {
        if (socket_wait(ftp->ctrl_fd, 0) <= 0) return -1; // socket is not ready for writing
 static int kftp_send_cmd(knetFile *ftp, const char *cmd, int is_get)
 {
        if (socket_wait(ftp->ctrl_fd, 0) <= 0) return -1; // socket is not ready for writing

-       write(ftp->ctrl_fd, cmd, strlen(cmd));
+       netwrite(ftp->ctrl_fd, cmd, strlen(cmd));

        return is_get? kftp_get_response(ftp) : 0;
 }
 
        return is_get? kftp_get_response(ftp) : 0;
 }
 


@@ -71,65 +226,39 @@ static int kftp_pasv_prep(knetFile *ftp)
        return 0;
 }
 
        return 0;
 }
 

+

 static int kftp_pasv_connect(knetFile *ftp)
 {
 static int kftp_pasv_connect(knetFile *ftp)
 {

-#define __err_pasv_connect(func) do { perror(func); freeaddrinfo(res); return -1; } while (0)
-
-       struct addrinfo hints, *res;
-       struct linger lng = { 0, 0 };
-       int on = 1;

        char host[80], port[10];
        char host[80], port[10];

-

        if (ftp->pasv_port == 0) {
                fprintf(stderr, "[kftp_pasv_connect] kftp_pasv_prep() is not called before hand.\n");
                return -1;
        }
        if (ftp->pasv_port == 0) {
                fprintf(stderr, "[kftp_pasv_connect] kftp_pasv_prep() is not called before hand.\n");
                return -1;
        }

-       memset(&hints, 0, sizeof(struct addrinfo));
-       hints.ai_family = AF_UNSPEC;
-       hints.ai_socktype = SOCK_STREAM;

        sprintf(host, "%d.%d.%d.%d", ftp->pasv_ip[0], ftp->pasv_ip[1], ftp->pasv_ip[2], ftp->pasv_ip[3]);
        sprintf(port, "%d", ftp->pasv_port);
        sprintf(host, "%d.%d.%d.%d", ftp->pasv_ip[0], ftp->pasv_ip[1], ftp->pasv_ip[2], ftp->pasv_ip[3]);
        sprintf(port, "%d", ftp->pasv_port);

-       if (getaddrinfo(host, port, &hints, &res) != 0) { perror("getaddrinfo"); return -1; }
-       if ((ftp->fd = socket(res->ai_family, res->ai_socktype, res->ai_protocol)) == -1) __err_pasv_connect("socket");
-       if (setsockopt(ftp->fd, SOL_SOCKET, SO_REUSEADDR, &on, sizeof(on)) == -1) __err_pasv_connect("setsockopt");
-       if (setsockopt(ftp->fd, SOL_SOCKET, SO_LINGER, &lng, sizeof(lng)) == -1) __err_pasv_connect("setsockopt");
-       if (connect(ftp->fd, res->ai_addr, res->ai_addrlen) != 0) __err_pasv_connect("connect");
-       freeaddrinfo(res);
+       ftp->fd = socket_connect(host, port);
+       if (ftp->fd == -1) return -1;

        return 0;
 }
 
 int kftp_connect(knetFile *ftp)
 {
        return 0;
 }
 
 int kftp_connect(knetFile *ftp)
 {

-#define __err_connect(func) do { perror(func); return -1; } while (0)
-
-       int on = 1;
-       { // open socket
-               struct addrinfo hints, *res;
-               memset(&hints, 0, sizeof(struct addrinfo));
-               hints.ai_family = AF_UNSPEC;
-               hints.ai_socktype = SOCK_STREAM;
-               if (getaddrinfo(ftp->host, "21", &hints, &res) != 0) __err_connect("getaddrinfo");
-               if ((ftp->ctrl_fd = socket(res->ai_family, res->ai_socktype, res->ai_protocol)) == -1) __err_connect("socket");
-               if (setsockopt(ftp->ctrl_fd, SOL_SOCKET, SO_REUSEADDR, &on, sizeof(on)) == -1) __err_connect("setsockopt");
-               if (connect(ftp->ctrl_fd, res->ai_addr, res->ai_addrlen) != 0) __err_connect("connect");
-               freeaddrinfo(res);
-               kftp_get_response(ftp);
-       }
-       { // login
-               kftp_send_cmd(ftp, "USER anonymous\r\n", 1);
-               kftp_send_cmd(ftp, "PASS kftp@\r\n", 1);
-               kftp_send_cmd(ftp, "TYPE I\r\n", 1);
-       }
+       ftp->ctrl_fd = socket_connect(ftp->host, ftp->port);
+       if (ftp->ctrl_fd == -1) return -1;
+       kftp_get_response(ftp);
+       kftp_send_cmd(ftp, "USER anonymous\r\n", 1);
+       kftp_send_cmd(ftp, "PASS kftp@\r\n", 1);
+       kftp_send_cmd(ftp, "TYPE I\r\n", 1);

        return 0;
 }
 
 int kftp_reconnect(knetFile *ftp)
 {
        return 0;
 }
 
 int kftp_reconnect(knetFile *ftp)
 {

-       if (ftp->ctrl_fd >= 0) {
-               close(ftp->ctrl_fd);
+       if (ftp->ctrl_fd != -1) {
+               netclose(ftp->ctrl_fd);

                ftp->ctrl_fd = -1;
        }
                ftp->ctrl_fd = -1;
        }

-       close(ftp->fd);
+       netclose(ftp->fd);

        return kftp_connect(ftp);
 }
 
        return kftp_connect(ftp);
 }
 


@@ -146,6 +275,9 @@ knetFile *kftp_parse_url(const char *fn, const char *mode)
        fp = calloc(1, sizeof(knetFile));
        fp->type = KNF_TYPE_FTP;
        fp->fd = -1;
        fp = calloc(1, sizeof(knetFile));
        fp->type = KNF_TYPE_FTP;
        fp->fd = -1;

+       /* the Linux/Mac version of socket_connect() also recognizes a port
+        * like "ftp", but the Windows version does not. */
+       fp->port = strdup("21");

        fp->host = calloc(l + 1, 1);
        if (strchr(mode, 'c')) fp->no_reconnect = 1;
        strncpy(fp->host, fn + 6, l);
        fp->host = calloc(l + 1, 1);
        if (strchr(mode, 'c')) fp->no_reconnect = 1;
        strncpy(fp->host, fn + 6, l);


@@ -158,14 +290,20 @@ knetFile *kftp_parse_url(const char *fn, const char *mode)
 int kftp_connect_file(knetFile *fp)
 {
        int ret;
 int kftp_connect_file(knetFile *fp)
 {
        int ret;

-       if (fp->fd >= 0) {
-               close(fp->fd);
+       if (fp->fd != -1) {
+               netclose(fp->fd);

                if (fp->no_reconnect) kftp_get_response(fp);
        }
        kftp_pasv_prep(fp);
        if (fp->offset) {
                char tmp[32];
                if (fp->no_reconnect) kftp_get_response(fp);
        }
        kftp_pasv_prep(fp);
        if (fp->offset) {
                char tmp[32];

+#ifndef _WIN32

                sprintf(tmp, "REST %lld\r\n", (long long)fp->offset);
                sprintf(tmp, "REST %lld\r\n", (long long)fp->offset);

+#else
+               strcpy(tmp, "REST ");
+               uint64tostr(tmp + 5, fp->offset);
+               strcat(tmp, "\r\n");
+#endif

                kftp_send_cmd(fp, tmp, 1);
        }
        kftp_send_cmd(fp, fp->retr, 0);
                kftp_send_cmd(fp, tmp, 1);
        }
        kftp_send_cmd(fp, fp->retr, 0);


@@ -173,7 +311,7 @@ int kftp_connect_file(knetFile *fp)
        ret = kftp_get_response(fp);
        if (ret != 150) {
                fprintf(stderr, "[kftp_connect_file] %s\n", fp->response);
        ret = kftp_get_response(fp);
        if (ret != 150) {
                fprintf(stderr, "[kftp_connect_file] %s\n", fp->response);

-               close(fp->fd);
+               netclose(fp->fd);

                fp->fd = -1;
                return -1;
        }
                fp->fd = -1;
                return -1;
        }


@@ -181,6 +319,92 @@ int kftp_connect_file(knetFile *fp)
        return 0;
 }
 
        return 0;
 }
 

+/**************************
+ * HTTP specific routines *
+ **************************/
+
+knetFile *khttp_parse_url(const char *fn, const char *mode)
+{
+       knetFile *fp;
+       char *p, *proxy, *q;
+       int l;
+       if (strstr(fn, "http://") != fn) return 0;
+       // set ->http_host
+       for (p = (char*)fn + 7; *p && *p != '/'; ++p);
+       l = p - fn - 7;
+       fp = calloc(1, sizeof(knetFile));
+       fp->http_host = calloc(l + 1, 1);
+       strncpy(fp->http_host, fn + 7, l);
+       fp->http_host[l] = 0;
+       for (q = fp->http_host; *q && *q != ':'; ++q);
+       if (*q == ':') *q++ = 0;
+       // get http_proxy
+       proxy = getenv("http_proxy");
+       // set ->host, ->port and ->path
+       if (proxy == 0) {
+               fp->host = strdup(fp->http_host); // when there is no proxy, server name is identical to http_host name.
+               fp->port = strdup(*q? q : "80");
+               fp->path = strdup(*p? p : "/");
+       } else {
+               fp->host = (strstr(proxy, "http://") == proxy)? strdup(proxy + 7) : strdup(proxy);
+               for (q = fp->host; *q && *q != ':'; ++q);
+               if (*q == ':') *q++ = 0; 
+               fp->port = strdup(*q? q : "80");
+               fp->path = strdup(fn);
+       }
+       fp->type = KNF_TYPE_HTTP;
+       fp->ctrl_fd = fp->fd = -1;
+       fp->seek_offset = -1;
+       return fp;
+}
+
+int khttp_connect_file(knetFile *fp)
+{
+       int ret, l = 0;
+       char *buf, *p;
+       if (fp->fd != -1) netclose(fp->fd);
+       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, "\r\n");
+       netwrite(fp->fd, buf, l);
+       l = 0;
+       while (netread(fp->fd, buf + l, 1)) { // read HTTP header; FIXME: bad efficiency
+               if (buf[l] == '\n' && l >= 3)
+                       if (strncmp(buf + l - 3, "\r\n\r\n", 4) == 0) break;
+               ++l;
+       }
+       buf[l] = 0;
+       if (l < 14) { // prematured header
+               netclose(fp->fd);
+               fp->fd = -1;
+               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
+               off_t rest = fp->offset;
+               while (rest) {
+                       off_t l = rest < 0x10000? rest : 0x10000;
+                       rest -= my_netread(fp->fd, buf, l);
+               }
+       } else if (ret != 206 && ret != 200) {
+               free(buf);
+               fprintf(stderr, "[khttp_connect_file] fail to open file (HTTP code: %d).\n", ret);
+               netclose(fp->fd);
+               fp->fd = -1;
+               return -1;
+       }
+       free(buf);
+       fp->is_ready = 1;
+       return 0;
+}
+
+/********************
+ * Generic routines *
+ ********************/
+

 knetFile *knet_open(const char *fn, const char *mode)
 {
        knetFile *fp = 0;
 knetFile *knet_open(const char *fn, const char *mode)
 {
        knetFile *fp = 0;


@@ -196,12 +420,19 @@ knetFile *knet_open(const char *fn, const char *mode)
                        return 0;
                }
                kftp_connect_file(fp);
                        return 0;
                }
                kftp_connect_file(fp);

-               if (fp->fd < 0) {
-                       knet_close(fp);
-                       return 0;
-               }
-       } else {
+       } else if (strstr(fn, "http://") == fn) {
+               fp = khttp_parse_url(fn, mode);
+               if (fp == 0) return 0;
+               khttp_connect_file(fp);
+       } else { // local file
+#ifdef _WIN32
+               /* In windows, O_BINARY is necessary. In Linux/Mac, O_BINARY may
+                * be undefined on some systems, although it is defined on my
+                * Mac and the Linux I have tested on. */
+               int fd = open(fn, O_RDONLY | O_BINARY);
+#else          

                int fd = open(fn, O_RDONLY);
                int fd = open(fn, O_RDONLY);

+#endif

                if (fd == -1) {
                        perror("open");
                        return 0;
                if (fd == -1) {
                        perror("open");
                        return 0;


@@ -209,6 +440,11 @@ knetFile *knet_open(const char *fn, const char *mode)
                fp = (knetFile*)calloc(1, sizeof(knetFile));
                fp->type = KNF_TYPE_LOCAL;
                fp->fd = fd;
                fp = (knetFile*)calloc(1, sizeof(knetFile));
                fp->type = KNF_TYPE_LOCAL;
                fp->fd = fd;

+               fp->ctrl_fd = -1;
+       }
+       if (fp && fp->fd == -1) {
+               knet_close(fp);
+               return 0;

        }
        return fp;
 }
        }
        return fp;
 }


@@ -224,46 +460,44 @@ knetFile *knet_dopen(int fd, const char *mode)
 off_t knet_read(knetFile *fp, void *buf, off_t len)
 {
        off_t l = 0;
 off_t knet_read(knetFile *fp, void *buf, off_t len)
 {
        off_t l = 0;

-       if (fp->fd < 0) return 0;
-       if (fp->type == KNF_TYPE_LOCAL) {
-               off_t rest = len, curr;
-               while (rest) {
-                       curr = read(fp->fd, buf + l, rest);
-                       if (curr == 0) break;
-                       l += curr; rest -= curr;
-               }
-               fp->offset += l;
-       } else {
-               off_t rest = len, curr;
+       if (fp->fd == -1) return 0;
+       if (fp->type == KNF_TYPE_FTP) {

                if (fp->is_ready == 0) {
                        if (!fp->no_reconnect) kftp_reconnect(fp);
                        kftp_connect_file(fp);
                if (fp->is_ready == 0) {
                        if (!fp->no_reconnect) kftp_reconnect(fp);
                        kftp_connect_file(fp);

-                       fp->is_ready = 1;

                }
                }

+       } else if (fp->type == KNF_TYPE_HTTP) {
+               if (fp->is_ready == 0)
+                       khttp_connect_file(fp);
+       }
+       if (fp->type == KNF_TYPE_LOCAL) { // on Windows, the following block is necessary; not on UNIX
+               off_t rest = len, curr;

                while (rest) {
                while (rest) {

-                       if (socket_wait(fp->fd, 1) <= 0) break; // socket is not ready for reading

                        curr = read(fp->fd, buf + l, rest);
                        curr = read(fp->fd, buf + l, rest);

-                       if (curr == 0) break; // FIXME: end of file or bad network? I do not know...
+                       if (curr == 0) break;

                        l += curr; rest -= curr;
                }
                        l += curr; rest -= curr;
                }

-               fp->offset += l;
-       }
+       } else l = my_netread(fp->fd, buf, len);
+       fp->offset += l;

        return l;
 }
 
 int knet_seek(knetFile *fp, off_t off, int whence)
 {
        return l;
 }
 
 int knet_seek(knetFile *fp, off_t off, int whence)
 {

+       if (whence == SEEK_SET && off == fp->offset) return 0;

        if (fp->type == KNF_TYPE_LOCAL) {
        if (fp->type == KNF_TYPE_LOCAL) {

-               if (lseek(fp->fd, off, whence) == -1) {
+               /* Be aware that lseek() returns the offset after seeking,
+                * while fseek() returns zero on success. */
+               off_t offset = lseek(fp->fd, off, whence);
+               if (offset == -1) {

                        perror("lseek");
                        return -1;
                }
                        perror("lseek");
                        return -1;
                }

-               fp->offset = off;
+               fp->offset = offset;

                return 0;
                return 0;

-       }
-       if (fp->type == KNF_TYPE_FTP) {
+       } 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
                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. Offset is unchanged.\n");
+                       fprintf(stderr, "[knet_seek] only SEEK_SET is supported for FTP/HTTP. Offset is unchanged.\n");

                        return -1;
                }
                fp->offset = off;
                        return -1;
                }
                fp->offset = off;


@@ -276,9 +510,16 @@ int knet_seek(knetFile *fp, off_t off, int whence)
 int knet_close(knetFile *fp)
 {
        if (fp == 0) return 0;
 int knet_close(knetFile *fp)
 {
        if (fp == 0) return 0;

-       if (fp->ctrl_fd >= 0) close(fp->ctrl_fd);
-       if (fp->fd >= 0) close(fp->fd);
-       free(fp->response); free(fp->retr); free(fp->host);
+       if (fp->ctrl_fd != -1) netclose(fp->ctrl_fd); // FTP specific
+       if (fp->fd != -1) {
+               /* On Linux/Mac, netclose() is an alias of close(), but on
+                * Windows, it is an alias of closesocket(). */
+               if (fp->type == KNF_TYPE_LOCAL) close(fp->fd);
+               else netclose(fp->fd);
+       }
+       free(fp->host); free(fp->port);
+       free(fp->response); free(fp->retr); // FTP specific
+       free(fp->path); free(fp->http_host); // HTTP specific

        free(fp);
        return 0;
 }
        free(fp);
        return 0;
 }


@@ -286,15 +527,40 @@ int knet_close(knetFile *fp)
 #ifdef KNETFILE_MAIN
 int main(void)
 {
 #ifdef KNETFILE_MAIN
 int main(void)
 {

-       char buf[256];
+       char *buf;

        knetFile *fp;
        knetFile *fp;

-//     fp = knet_open("ftp://ftp.ncbi.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom6.SLX.SRP000032.2009_06.bam", "r"); knet_seek(fp, 2500000000ll, SEEK_SET);
-       fp = knet_open("ftp://ftp.sanger.ac.uk/pub4/treefam/tmp/index.shtml", "r"); knet_seek(fp, 2000, SEEK_SET);
-//     fp = knet_open("knetfile.c", "r"); knet_seek(fp, 2000, SEEK_SET);
-       knet_read(fp, buf, 255);
-       buf[255] = 0;
-       printf("%s\n", buf);
+       int type = 4, l;
+#ifdef _WIN32
+       knet_win32_init();
+#endif
+       buf = calloc(0x100000, 1);
+       if (type == 0) {
+               fp = knet_open("knetfile.c", "r");
+               knet_seek(fp, 1000, SEEK_SET);
+       } else if (type == 1) { // NCBI FTP, large file
+               fp = knet_open("ftp://ftp.ncbi.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom6.SLX.SRP000032.2009_06.bam", "r");
+               knet_seek(fp, 2500000000ll, SEEK_SET);
+               l = knet_read(fp, buf, 255);
+       } else if (type == 2) {
+               fp = knet_open("ftp://ftp.sanger.ac.uk/pub4/treefam/tmp/index.shtml", "r");
+               knet_seek(fp, 1000, SEEK_SET);
+       } else if (type == 3) {
+               fp = knet_open("http://www.sanger.ac.uk/Users/lh3/index.shtml", "r");
+               knet_seek(fp, 1000, SEEK_SET);
+       } else if (type == 4) {
+               fp = knet_open("http://www.sanger.ac.uk/Users/lh3/ex1.bam", "r");
+               knet_read(fp, buf, 10000);
+               knet_seek(fp, 20000, SEEK_SET);
+               knet_seek(fp, 10000, SEEK_SET);
+               l = knet_read(fp, buf+10000, 10000000) + 10000;
+       }
+       if (type != 4 && type != 1) {
+               knet_read(fp, buf, 255);
+               buf[255] = 0;
+               printf("%s\n", buf);
+       } else write(fileno(stdout), buf, l);

        knet_close(fp);
        knet_close(fp);

+       free(buf);

        return 0;
 }
 #endif
        return 0;
 }
 #endif

diff --git a/knetfile.h b/knetfile.h
index bf45f3df809b8f3ebd54b17c711dd188ab3dbb36..9021b937ac211c1459ef1afec3c20b49ce6007a1 100644 (file)
--- a/knetfile.h
+++ b/knetfile.h
@@ -4,6 +4,17 @@
 #include <stdint.h>
 #include <fcntl.h>
 
 #include <stdint.h>
 #include <fcntl.h>
 

+#ifndef _WIN32
+#define netread(fd, ptr, len) read(fd, ptr, len)
+#define netwrite(fd, ptr, len) write(fd, ptr, len)
+#define netclose(fd) close(fd)
+#else
+#include <winsock.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)
+#endif
+

 // FIXME: currently I/O is unbuffered
 
 #define KNF_TYPE_LOCAL 1
 // FIXME: currently I/O is unbuffered
 
 #define KNF_TYPE_LOCAL 1


@@ -13,12 +24,15 @@
 typedef struct knetFile_s {
        int type, fd;
        int64_t offset;
 typedef struct knetFile_s {
        int type, fd;
        int64_t offset;

-       char *host;
+       char *host, *port;

 
        // the following are for FTP only
        int ctrl_fd, pasv_ip[4], pasv_port, max_response, no_reconnect, is_ready;
        char *response, *retr;
        int64_t seek_offset; // for lazy seek
 
        // the following are for FTP only
        int ctrl_fd, pasv_ip[4], pasv_port, max_response, no_reconnect, is_ready;
        char *response, *retr;
        int64_t seek_offset; // for lazy seek

+
+       // the following are for HTTP only
+       char *path, *http_host;

 } knetFile;
 
 #define knet_tell(fp) ((fp)->offset)
 } knetFile;
 
 #define knet_tell(fp) ((fp)->offset)


@@ -28,6 +42,11 @@ typedef struct knetFile_s {
 extern "C" {
 #endif
 
 extern "C" {
 #endif
 

+#ifdef _WIN32
+       int knet_win32_init();
+       void knet_win32_destroy();
+#endif
+

        knetFile *knet_open(const char *fn, const char *mode);
 
        /* 
        knetFile *knet_open(const char *fn, const char *mode);
 
        /* 

diff --git a/kseq.h b/kseq.h
index bbe012552ff3b2ac047ddd10655e80827bbb4008..82face095919a3991b1f927bb5422ffd95f102a4 100644 (file)
--- a/kseq.h
+++ b/kseq.h
@@ -25,6 +25,10 @@
 
 /* Contact: Heng Li <lh3@sanger.ac.uk> */
 
 
 /* Contact: Heng Li <lh3@sanger.ac.uk> */
 

+/*
+  2009-07-16 (lh3): in kstream_t, change "char*" to "unsigned char*"
+ */
+

 /* Last Modified: 12APR2009 */
 
 #ifndef AC_KSEQ_H
 /* Last Modified: 12APR2009 */
 
 #ifndef AC_KSEQ_H


@@ -40,7 +44,7 @@
 
 #define __KS_TYPE(type_t)                                              \
        typedef struct __kstream_t {                            \
 
 #define __KS_TYPE(type_t)                                              \
        typedef struct __kstream_t {                            \

-               char *buf;                                                              \
+               unsigned char *buf;                                             \

                int begin, end, is_eof;                                 \
                type_t f;                                                               \
        } kstream_t;
                int begin, end, is_eof;                                 \
                type_t f;                                                               \
        } kstream_t;


@@ -53,7 +57,7 @@
        {                                                                                                                               \
                kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t));       \
                ks->f = f;                                                                                                      \
        {                                                                                                                               \
                kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t));       \
                ks->f = f;                                                                                                      \

-               ks->buf = (char*)malloc(__bufsize);                                                     \
+               ks->buf = malloc(__bufsize);                                                            \

                return ks;                                                                                                      \
        }                                                                                                                               \
        static inline void ks_destroy(kstream_t *ks)                                    \
                return ks;                                                                                                      \
        }                                                                                                                               \
        static inline void ks_destroy(kstream_t *ks)                                    \

diff --git a/kstring.c b/kstring.c
index dc20faeac9fef21c88813c34f9f42ac6c69d449d..e0203fa3e877604f53291b39ff6a2f273e6dff8f 100644 (file)
--- a/kstring.c
+++ b/kstring.c
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <ctype.h>
 #include <string.h>
 #include <stdio.h>
 #include <ctype.h>
 #include <string.h>

+#include <stdint.h>

 #include "kstring.h"
 
 int ksprintf(kstring_t *s, const char *fmt, ...)
 #include "kstring.h"
 
 int ksprintf(kstring_t *s, const char *fmt, ...)


@@ -61,6 +62,78 @@ int ksplit_core(char *s, int delimiter, int *_max, int **_offsets)
        return n;
 }
 
        return n;
 }
 

+/**********************
+ * Boyer-Moore search *
+ **********************/
+
+// reference: http://www-igm.univ-mlv.fr/~lecroq/string/node14.html
+int *ksBM_prep(const uint8_t *pat, int m)
+{
+       int i, *suff, *prep, *bmGs, *bmBc;
+       prep = calloc(m + 256, 1);
+       bmGs = prep; bmBc = prep + m;
+       { // preBmBc()
+               for (i = 0; i < 256; ++i) bmBc[i] = m;
+               for (i = 0; i < m - 1; ++i) bmBc[pat[i]] = m - i - 1;
+       }
+       suff = calloc(m, sizeof(int));
+       { // suffixes()
+               int f = 0, g;
+               suff[m - 1] = m;
+               g = m - 1;
+               for (i = m - 2; i >= 0; --i) {
+                       if (i > g && suff[i + m - 1 - f] < i - g)
+                               suff[i] = suff[i + m - 1 - f];
+                       else {
+                               if (i < g) g = i;
+                               f = i;
+                               while (g >= 0 && pat[g] == pat[g + m - 1 - f]) --g;
+                               suff[i] = f - g;
+                       }
+               }
+       }
+       { // preBmGs()
+               int j = 0;
+               for (i = 0; i < m; ++i) bmGs[i] = m;
+               for (i = m - 1; i >= 0; --i)
+                       if (suff[i] == i + 1)
+                               for (; j < m - 1 - i; ++j)
+                                       if (bmGs[j] == m)
+                                               bmGs[j] = m - 1 - i;
+               for (i = 0; i <= m - 2; ++i)
+                       bmGs[m - 1 - suff[i]] = m - 1 - i;
+       }
+       free(suff);
+       return prep;
+}
+
+int *ksBM_search(const uint8_t *str, int n, const uint8_t *pat, int m, int *_prep, int *n_matches)
+{
+       int i, j, *prep, *bmGs, *bmBc;
+       int *matches = 0, mm = 0, nm = 0;
+       prep = _prep? _prep : ksBM_prep(pat, m);
+       bmGs = prep; bmBc = prep + m;
+       j = 0;
+       while (j <= n - m) {
+               for (i = m - 1; i >= 0 && pat[i] == str[i+j]; --i);
+               if (i < 0) {
+                       if (nm == mm) {
+                               mm = mm? mm<<1 : 1;
+                               matches = realloc(matches, mm * sizeof(int));
+                       }
+                       matches[nm++] = j;
+                       j += bmGs[0];
+               } else {
+                       int max = bmBc[str[i+j]] - m + 1 + i;
+                       if (max < bmGs[i]) max = bmGs[i];
+                       j += max;
+               }
+       }
+       *n_matches = nm;
+       if (_prep == 0) free(prep);
+       return matches;
+}
+

 #ifdef KSTRING_MAIN
 #include <stdio.h>
 int main()
 #ifdef KSTRING_MAIN
 #include <stdio.h>
 int main()


@@ -76,6 +149,17 @@ int main()
        for (i = 0; i < n; ++i)
                printf("field[%d] = '%s'\n", i, s->s + fields[i]);
        free(s);
        for (i = 0; i < n; ++i)
                printf("field[%d] = '%s'\n", i, s->s + fields[i]);
        free(s);

+
+       {
+               static char *str = "abcdefgcdg";
+               static char *pat = "cd";
+               int n, *matches;
+               matches = ksBM_search(str, strlen(str), pat, strlen(pat), 0, &n);
+               printf("%d: \n", n);
+               for (i = 0; i < n; ++i)
+                       printf("- %d\n", matches[i]);
+               free(matches);
+       }

        return 0;
 }
 #endif
        return 0;
 }
 #endif

diff --git a/kstring.h b/kstring.h
index 221ade2472655ba842477b1249feb7ca7209631c..f4e5a99df5b03364166338e55799edbc41b1cae4 100644 (file)
--- a/kstring.h
+++ b/kstring.h
@@ -3,6 +3,7 @@
 
 #include <stdlib.h>
 #include <string.h>
 
 #include <stdlib.h>
 #include <string.h>

+#include <stdint.h>

 
 #ifndef kroundup32
 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
 
 #ifndef kroundup32
 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))


@@ -19,6 +20,14 @@ typedef struct __kstring_t {
 int ksprintf(kstring_t *s, const char *fmt, ...);
 int ksplit_core(char *s, int delimiter, int *_max, int **_offsets);
 
 int ksprintf(kstring_t *s, const char *fmt, ...);
 int ksplit_core(char *s, int delimiter, int *_max, int **_offsets);
 

+// calculate the auxiliary array, allocated by calloc()
+int *ksBM_prep(const uint8_t *pat, int m);
+
+/* Search pat in str and returned the list of matches. The size of the
+ * list is returned as n_matches. _prep is the array returned by
+ * ksBM_prep(). If it is a NULL pointer, ksBM_prep() will be called. */
+int *ksBM_search(const uint8_t *str, int n, const uint8_t *pat, int m, int *_prep, int *n_matches);
+

 static inline int kputsn(const char *p, int l, kstring_t *s)
 {
        if (s->l + l + 1 >= s->m) {
 static inline int kputsn(const char *p, int l, kstring_t *s)
 {
        if (s->l + l + 1 >= s->m) {

diff --git a/misc/maq2sam.c b/misc/maq2sam.c
index 758a698ed59aba300720599cda96c2759317887f..2bfbe2a3b328e276b52f1a443dbc6500d81122e8 100644 (file)
--- a/misc/maq2sam.c
+++ b/misc/maq2sam.c
@@ -5,7 +5,7 @@
 #include <stdlib.h>
 #include <assert.h>
 
 #include <stdlib.h>
 #include <assert.h>
 

-#define PACKAGE_VERSION "0.1.2 (20090521)"
+#define PACKAGE_VERSION "r439"

 
 //#define MAQ_LONGREADS
 
 
 //#define MAQ_LONGREADS
 


@@ -111,9 +111,9 @@ void maq2tam_core(gzFile fp, const char *rg)
                                else if (m1->flag&(PAIRFLAG_RF|PAIRFLAG_RR)) c = 1;
                                else c = m1->pos&1;
                        }
                                else if (m1->flag&(PAIRFLAG_RF|PAIRFLAG_RR)) c = 1;
                                else c = m1->pos&1;
                        }

-                       flag |= c;
+                       if (c) flag |= 0x20;

                }
                }

-               if (flag) {
+               if (m1->flag) {

                        int l = strlen(m1->name);
                        if (m1->name[l-2] == '/') {
                                flag |= (m1->name[l-1] == '1')? 0x40 : 0x80;
                        int l = strlen(m1->name);
                        if (m1->name[l-2] == '/') {
                                flag |= (m1->name[l-1] == '1')? 0x40 : 0x80;

diff --git a/misc/md5.c b/misc/md5.c
index ccead0e6b8937f28ea778b4d80766db83f7839b6..55ae181c82d2a0e94c1a443799a7a72827c6df33 100644 (file)
--- a/misc/md5.c
+++ b/misc/md5.c
@@ -1,271 +1,260 @@
 /*
 /*

- **********************************************************************
- ** md5.c                                                            **
- ** RSA Data Security, Inc. MD5 Message Digest Algorithm             **
- ** Created: 2/17/90 RLR                                             **
- ** Revised: 1/91 SRD,AJ,BSK,JT Reference C Version                  **
- **********************************************************************
+ * This code implements the MD5 message-digest algorithm.
+ * The algorithm is due to Ron Rivest. This code was
+ * written by Colin Plumb in 1993, no copyright is claimed.
+ * This code is in the public domain; do with it what you wish.
+ *
+ * Equivalent code is available from RSA Data Security, Inc.
+ * This code has been tested against that, and is equivalent,
+ * except that you don't need to include two pages of legalese
+ * with every copy.
+ *
+ * To compute the message digest of a chunk of bytes, declare an
+ * MD5Context structure, pass it to MD5Init, call MD5Update as
+ * needed on buffers full of bytes, and then call MD5Final, which
+ * will fill a supplied 16-byte array with the digest.

  */
 
  */
 

+/* Brutally hacked by John Walker back from ANSI C to K&R (no
+   prototypes) to maintain the tradition that Netfone will compile
+   with Sun's original "cc". */
+
+#include <string.h>
+#include "md5.h"
+
+#ifndef HIGHFIRST
+#define byteReverse(buf, len)  /* Nothing */
+#else

 /*
 /*

- **********************************************************************
- ** Copyright (C) 1990, RSA Data Security, Inc. All rights reserved. **
- **                                                                  **
- ** License to copy and use this software is granted provided that   **
- ** it is identified as the "RSA Data Security, Inc. MD5 Message     **
- ** Digest Algorithm" in all material mentioning or referencing this **
- ** software or this function.                                       **
- **                                                                  **
- ** License is also granted to make and use derivative works         **
- ** provided that such works are identified as "derived from the RSA **
- ** Data Security, Inc. MD5 Message Digest Algorithm" in all         **
- ** material mentioning or referencing the derived work.             **
- **                                                                  **
- ** RSA Data Security, Inc. makes no representations concerning      **
- ** either the merchantability of this software or the suitability   **
- ** of this software for any particular purpose.  It is provided "as **
- ** is" without express or implied warranty of any kind.             **
- **                                                                  **
- ** These notices must be retained in any copies of any part of this **
- ** documentation and/or software.                                   **
- **********************************************************************
+ * Note: this code is harmless on little-endian machines.

  */
  */

+void byteReverse(buf, longs)
+    unsigned char *buf; unsigned longs;
+{
+    uint32_t t;
+    do {
+       t = (uint32_t) ((unsigned) buf[3] << 8 | buf[2]) << 16 |
+           ((unsigned) buf[1] << 8 | buf[0]);
+       *(uint32_t *) buf = t;
+       buf += 4;
+    } while (--longs);
+}
+#endif
+
+void MD5Transform(uint32_t buf[4], uint32_t in[16]);

 
 

-#include "md5.h"

 
 

-/* forward declaration */
-static void Transform ();
-
-static unsigned char PADDING[64] = {
-  0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
-};
-
-/* F, G and H are basic MD5 functions: selection, majority, parity */
-#define F(x, y, z) (((x) & (y)) | ((~x) & (z)))
-#define G(x, y, z) (((x) & (z)) | ((y) & (~z)))
-#define H(x, y, z) ((x) ^ (y) ^ (z))
-#define I(x, y, z) ((y) ^ ((x) | (~z))) 
-
-/* ROTATE_LEFT rotates x left n bits */
-#define ROTATE_LEFT(x, n) (((x) << (n)) | ((x) >> (32-(n))))
-
-/* FF, GG, HH, and II transformations for rounds 1, 2, 3, and 4 */
-/* Rotation is separate from addition to prevent recomputation */
-#define FF(a, b, c, d, x, s, ac) \
-  {(a) += F ((b), (c), (d)) + (x) + (UINT4)(ac); \
-   (a) = ROTATE_LEFT ((a), (s)); \
-   (a) += (b); \
-  }
-#define GG(a, b, c, d, x, s, ac) \
-  {(a) += G ((b), (c), (d)) + (x) + (UINT4)(ac); \
-   (a) = ROTATE_LEFT ((a), (s)); \
-   (a) += (b); \
-  }
-#define HH(a, b, c, d, x, s, ac) \
-  {(a) += H ((b), (c), (d)) + (x) + (UINT4)(ac); \
-   (a) = ROTATE_LEFT ((a), (s)); \
-   (a) += (b); \
-  }
-#define II(a, b, c, d, x, s, ac) \
-  {(a) += I ((b), (c), (d)) + (x) + (UINT4)(ac); \
-   (a) = ROTATE_LEFT ((a), (s)); \
-   (a) += (b); \
-  }
-
-void MD5Init (mdContext)
-MD5_CTX *mdContext;
+/*
+ * Start MD5 accumulation.  Set bit count to 0 and buffer to mysterious
+ * initialization constants.
+ */
+void MD5Init(ctx)
+    struct MD5Context *ctx;

 {
 {

-  mdContext->i[0] = mdContext->i[1] = (UINT4)0;
-
-  /* Load magic initialization constants.
-   */
-  mdContext->buf[0] = (UINT4)0x67452301;
-  mdContext->buf[1] = (UINT4)0xefcdab89;
-  mdContext->buf[2] = (UINT4)0x98badcfe;
-  mdContext->buf[3] = (UINT4)0x10325476;
+    ctx->buf[0] = 0x67452301;
+    ctx->buf[1] = 0xefcdab89;
+    ctx->buf[2] = 0x98badcfe;
+    ctx->buf[3] = 0x10325476;
+
+    ctx->bits[0] = 0;
+    ctx->bits[1] = 0;

 }
 
 }
 

-void MD5Update (mdContext, inBuf, inLen)
-MD5_CTX *mdContext;
-unsigned char *inBuf;
-unsigned int inLen;
+/*
+ * Update context to reflect the concatenation of another buffer full
+ * of bytes.
+ */
+void MD5Update(ctx, buf, len)
+    struct MD5Context *ctx; unsigned char *buf; unsigned len;

 {
 {

-  UINT4 in[16];
-  int mdi;
-  unsigned int i, ii;
-
-  /* compute number of bytes mod 64 */
-  mdi = (int)((mdContext->i[0] >> 3) & 0x3F);
-
-  /* update number of bits */
-  if ((mdContext->i[0] + ((UINT4)inLen << 3)) < mdContext->i[0])
-    mdContext->i[1]++;
-  mdContext->i[0] += ((UINT4)inLen << 3);
-  mdContext->i[1] += ((UINT4)inLen >> 29);
-
-  while (inLen--) {
-    /* add new character to buffer, increment mdi */
-    mdContext->in[mdi++] = *inBuf++;
-
-    /* transform if necessary */
-    if (mdi == 0x40) {
-      for (i = 0, ii = 0; i < 16; i++, ii += 4)
-        in[i] = (((UINT4)mdContext->in[ii+3]) << 24) |
-                (((UINT4)mdContext->in[ii+2]) << 16) |
-                (((UINT4)mdContext->in[ii+1]) << 8) |
-                ((UINT4)mdContext->in[ii]);
-      Transform (mdContext->buf, in);
-      mdi = 0;
+    uint32_t t;
+
+    /* Update bitcount */
+
+    t = ctx->bits[0];
+    if ((ctx->bits[0] = t + ((uint32_t) len << 3)) < t)
+       ctx->bits[1]++;         /* Carry from low to high */
+    ctx->bits[1] += len >> 29;
+
+    t = (t >> 3) & 0x3f;       /* Bytes already in shsInfo->data */
+
+    /* Handle any leading odd-sized chunks */
+
+    if (t) {
+       unsigned char *p = (unsigned char *) ctx->in + t;
+
+       t = 64 - t;
+       if (len < t) {
+           memcpy(p, buf, len);
+           return;
+       }
+       memcpy(p, buf, t);
+       byteReverse(ctx->in, 16);
+       MD5Transform(ctx->buf, (uint32_t *) ctx->in);
+       buf += t;
+       len -= t;
+    }
+    /* Process data in 64-byte chunks */
+
+    while (len >= 64) {
+       memcpy(ctx->in, buf, 64);
+       byteReverse(ctx->in, 16);
+       MD5Transform(ctx->buf, (uint32_t *) ctx->in);
+       buf += 64;
+       len -= 64;

     }
     }

-  }
+
+    /* Handle any remaining bytes of data. */
+
+    memcpy(ctx->in, buf, len);

 }
 
 }
 

-void MD5Final (mdContext)
-MD5_CTX *mdContext;
+/*
+ * Final wrapup - pad to 64-byte boundary with the bit pattern 
+ * 1 0* (64-bit count of bits processed, MSB-first)
+ */
+void MD5Final(digest, ctx)
+    unsigned char digest[16]; struct MD5Context *ctx;

 {
 {

-  UINT4 in[16];
-  int mdi;
-  unsigned int i, ii;
-  unsigned int padLen;
-
-  /* save number of bits */
-  in[14] = mdContext->i[0];
-  in[15] = mdContext->i[1];
-
-  /* compute number of bytes mod 64 */
-  mdi = (int)((mdContext->i[0] >> 3) & 0x3F);
-
-  /* pad out to 56 mod 64 */
-  padLen = (mdi < 56) ? (56 - mdi) : (120 - mdi);
-  MD5Update (mdContext, PADDING, padLen);
-
-  /* append length in bits and transform */
-  for (i = 0, ii = 0; i < 14; i++, ii += 4)
-    in[i] = (((UINT4)mdContext->in[ii+3]) << 24) |
-            (((UINT4)mdContext->in[ii+2]) << 16) |
-            (((UINT4)mdContext->in[ii+1]) << 8) |
-            ((UINT4)mdContext->in[ii]);
-  Transform (mdContext->buf, in);
-
-  /* store buffer in digest */
-  for (i = 0, ii = 0; i < 4; i++, ii += 4) {
-    mdContext->digest[ii] = (unsigned char)(mdContext->buf[i] & 0xFF);
-    mdContext->digest[ii+1] =
-      (unsigned char)((mdContext->buf[i] >> 8) & 0xFF);
-    mdContext->digest[ii+2] =
-      (unsigned char)((mdContext->buf[i] >> 16) & 0xFF);
-    mdContext->digest[ii+3] =
-      (unsigned char)((mdContext->buf[i] >> 24) & 0xFF);
-  }
+    unsigned count;
+    unsigned char *p;
+
+    /* Compute number of bytes mod 64 */
+    count = (ctx->bits[0] >> 3) & 0x3F;
+
+    /* Set the first char of padding to 0x80.  This is safe since there is
+       always at least one byte free */
+    p = ctx->in + count;
+    *p++ = 0x80;
+
+    /* Bytes of padding needed to make 64 bytes */
+    count = 64 - 1 - count;
+
+    /* Pad out to 56 mod 64 */
+    if (count < 8) {
+       /* Two lots of padding:  Pad the first block to 64 bytes */
+       memset(p, 0, count);
+       byteReverse(ctx->in, 16);
+       MD5Transform(ctx->buf, (uint32_t *) ctx->in);
+
+       /* Now fill the next block with 56 bytes */
+       memset(ctx->in, 0, 56);
+    } else {
+       /* Pad block to 56 bytes */
+       memset(p, 0, count - 8);
+    }
+    byteReverse(ctx->in, 14);
+
+    /* Append length in bits and transform */
+    ((uint32_t *) ctx->in)[14] = ctx->bits[0];
+    ((uint32_t *) ctx->in)[15] = ctx->bits[1];
+
+    MD5Transform(ctx->buf, (uint32_t *) ctx->in);
+    byteReverse((unsigned char *) ctx->buf, 4);
+    memcpy(digest, ctx->buf, 16);
+    memset(ctx, 0, sizeof(ctx));        /* In case it's sensitive */

 }
 
 }
 

-/* Basic MD5 step. Transform buf based on in.
+
+/* The four core functions - F1 is optimized somewhat */
+
+/* #define F1(x, y, z) (x & y | ~x & z) */
+#define F1(x, y, z) (z ^ (x & (y ^ z)))
+#define F2(x, y, z) F1(z, x, y)
+#define F3(x, y, z) (x ^ y ^ z)
+#define F4(x, y, z) (y ^ (x | ~z))
+
+/* This is the central step in the MD5 algorithm. */
+#define MD5STEP(f, w, x, y, z, data, s) \
+       ( w += f(x, y, z) + data,  w = w<<s | w>>(32-s),  w += x )
+
+/*
+ * The core of the MD5 algorithm, this alters an existing MD5 hash to
+ * reflect the addition of 16 longwords of new data.  MD5Update blocks
+ * the data and converts bytes into longwords for this routine.

  */
  */

-static void Transform (buf, in)
-UINT4 *buf;
-UINT4 *in;
+void MD5Transform(buf, in)
+    uint32_t buf[4]; uint32_t in[16];

 {
 {

-  UINT4 a = buf[0], b = buf[1], c = buf[2], d = buf[3];
-
-  /* Round 1 */
-#define S11 7
-#define S12 12
-#define S13 17
-#define S14 22
-  FF ( a, b, c, d, in[ 0], S11, 3614090360u); /* 1 */
-  FF ( d, a, b, c, in[ 1], S12, 3905402710u); /* 2 */
-  FF ( c, d, a, b, in[ 2], S13,  606105819u); /* 3 */
-  FF ( b, c, d, a, in[ 3], S14, 3250441966u); /* 4 */
-  FF ( a, b, c, d, in[ 4], S11, 4118548399u); /* 5 */
-  FF ( d, a, b, c, in[ 5], S12, 1200080426u); /* 6 */
-  FF ( c, d, a, b, in[ 6], S13, 2821735955u); /* 7 */
-  FF ( b, c, d, a, in[ 7], S14, 4249261313u); /* 8 */
-  FF ( a, b, c, d, in[ 8], S11, 1770035416u); /* 9 */
-  FF ( d, a, b, c, in[ 9], S12, 2336552879u); /* 10 */
-  FF ( c, d, a, b, in[10], S13, 4294925233u); /* 11 */
-  FF ( b, c, d, a, in[11], S14, 2304563134u); /* 12 */
-  FF ( a, b, c, d, in[12], S11, 1804603682u); /* 13 */
-  FF ( d, a, b, c, in[13], S12, 4254626195u); /* 14 */
-  FF ( c, d, a, b, in[14], S13, 2792965006u); /* 15 */
-  FF ( b, c, d, a, in[15], S14, 1236535329u); /* 16 */
-
-  /* Round 2 */
-#define S21 5
-#define S22 9
-#define S23 14
-#define S24 20
-  GG ( a, b, c, d, in[ 1], S21, 4129170786u); /* 17 */
-  GG ( d, a, b, c, in[ 6], S22, 3225465664u); /* 18 */
-  GG ( c, d, a, b, in[11], S23,  643717713u); /* 19 */
-  GG ( b, c, d, a, in[ 0], S24, 3921069994u); /* 20 */
-  GG ( a, b, c, d, in[ 5], S21, 3593408605u); /* 21 */
-  GG ( d, a, b, c, in[10], S22,   38016083u); /* 22 */
-  GG ( c, d, a, b, in[15], S23, 3634488961u); /* 23 */
-  GG ( b, c, d, a, in[ 4], S24, 3889429448u); /* 24 */
-  GG ( a, b, c, d, in[ 9], S21,  568446438u); /* 25 */
-  GG ( d, a, b, c, in[14], S22, 3275163606u); /* 26 */
-  GG ( c, d, a, b, in[ 3], S23, 4107603335u); /* 27 */
-  GG ( b, c, d, a, in[ 8], S24, 1163531501u); /* 28 */
-  GG ( a, b, c, d, in[13], S21, 2850285829u); /* 29 */
-  GG ( d, a, b, c, in[ 2], S22, 4243563512u); /* 30 */
-  GG ( c, d, a, b, in[ 7], S23, 1735328473u); /* 31 */
-  GG ( b, c, d, a, in[12], S24, 2368359562u); /* 32 */
-
-  /* Round 3 */
-#define S31 4
-#define S32 11
-#define S33 16
-#define S34 23
-  HH ( a, b, c, d, in[ 5], S31, 4294588738u); /* 33 */
-  HH ( d, a, b, c, in[ 8], S32, 2272392833u); /* 34 */
-  HH ( c, d, a, b, in[11], S33, 1839030562u); /* 35 */
-  HH ( b, c, d, a, in[14], S34, 4259657740u); /* 36 */
-  HH ( a, b, c, d, in[ 1], S31, 2763975236u); /* 37 */
-  HH ( d, a, b, c, in[ 4], S32, 1272893353u); /* 38 */
-  HH ( c, d, a, b, in[ 7], S33, 4139469664u); /* 39 */
-  HH ( b, c, d, a, in[10], S34, 3200236656u); /* 40 */
-  HH ( a, b, c, d, in[13], S31,  681279174u); /* 41 */
-  HH ( d, a, b, c, in[ 0], S32, 3936430074u); /* 42 */
-  HH ( c, d, a, b, in[ 3], S33, 3572445317u); /* 43 */
-  HH ( b, c, d, a, in[ 6], S34,   76029189u); /* 44 */
-  HH ( a, b, c, d, in[ 9], S31, 3654602809u); /* 45 */
-  HH ( d, a, b, c, in[12], S32, 3873151461u); /* 46 */
-  HH ( c, d, a, b, in[15], S33,  530742520u); /* 47 */
-  HH ( b, c, d, a, in[ 2], S34, 3299628645u); /* 48 */
-
-  /* Round 4 */
-#define S41 6
-#define S42 10
-#define S43 15
-#define S44 21
-  II ( a, b, c, d, in[ 0], S41, 4096336452u); /* 49 */
-  II ( d, a, b, c, in[ 7], S42, 1126891415u); /* 50 */
-  II ( c, d, a, b, in[14], S43, 2878612391u); /* 51 */
-  II ( b, c, d, a, in[ 5], S44, 4237533241u); /* 52 */
-  II ( a, b, c, d, in[12], S41, 1700485571u); /* 53 */
-  II ( d, a, b, c, in[ 3], S42, 2399980690u); /* 54 */
-  II ( c, d, a, b, in[10], S43, 4293915773u); /* 55 */
-  II ( b, c, d, a, in[ 1], S44, 2240044497u); /* 56 */
-  II ( a, b, c, d, in[ 8], S41, 1873313359u); /* 57 */
-  II ( d, a, b, c, in[15], S42, 4264355552u); /* 58 */
-  II ( c, d, a, b, in[ 6], S43, 2734768916u); /* 59 */
-  II ( b, c, d, a, in[13], S44, 1309151649u); /* 60 */
-  II ( a, b, c, d, in[ 4], S41, 4149444226u); /* 61 */
-  II ( d, a, b, c, in[11], S42, 3174756917u); /* 62 */
-  II ( c, d, a, b, in[ 2], S43,  718787259u); /* 63 */
-  II ( b, c, d, a, in[ 9], S44, 3951481745u); /* 64 */
-
-  buf[0] += a;
-  buf[1] += b;
-  buf[2] += c;
-  buf[3] += d;
+    register uint32_t a, b, c, d;
+
+    a = buf[0];
+    b = buf[1];
+    c = buf[2];
+    d = buf[3];
+
+    MD5STEP(F1, a, b, c, d, in[0] + 0xd76aa478, 7);
+    MD5STEP(F1, d, a, b, c, in[1] + 0xe8c7b756, 12);
+    MD5STEP(F1, c, d, a, b, in[2] + 0x242070db, 17);
+    MD5STEP(F1, b, c, d, a, in[3] + 0xc1bdceee, 22);
+    MD5STEP(F1, a, b, c, d, in[4] + 0xf57c0faf, 7);
+    MD5STEP(F1, d, a, b, c, in[5] + 0x4787c62a, 12);
+    MD5STEP(F1, c, d, a, b, in[6] + 0xa8304613, 17);
+    MD5STEP(F1, b, c, d, a, in[7] + 0xfd469501, 22);
+    MD5STEP(F1, a, b, c, d, in[8] + 0x698098d8, 7);
+    MD5STEP(F1, d, a, b, c, in[9] + 0x8b44f7af, 12);
+    MD5STEP(F1, c, d, a, b, in[10] + 0xffff5bb1, 17);
+    MD5STEP(F1, b, c, d, a, in[11] + 0x895cd7be, 22);
+    MD5STEP(F1, a, b, c, d, in[12] + 0x6b901122, 7);
+    MD5STEP(F1, d, a, b, c, in[13] + 0xfd987193, 12);
+    MD5STEP(F1, c, d, a, b, in[14] + 0xa679438e, 17);
+    MD5STEP(F1, b, c, d, a, in[15] + 0x49b40821, 22);
+
+    MD5STEP(F2, a, b, c, d, in[1] + 0xf61e2562, 5);
+    MD5STEP(F2, d, a, b, c, in[6] + 0xc040b340, 9);
+    MD5STEP(F2, c, d, a, b, in[11] + 0x265e5a51, 14);
+    MD5STEP(F2, b, c, d, a, in[0] + 0xe9b6c7aa, 20);
+    MD5STEP(F2, a, b, c, d, in[5] + 0xd62f105d, 5);
+    MD5STEP(F2, d, a, b, c, in[10] + 0x02441453, 9);
+    MD5STEP(F2, c, d, a, b, in[15] + 0xd8a1e681, 14);
+    MD5STEP(F2, b, c, d, a, in[4] + 0xe7d3fbc8, 20);
+    MD5STEP(F2, a, b, c, d, in[9] + 0x21e1cde6, 5);
+    MD5STEP(F2, d, a, b, c, in[14] + 0xc33707d6, 9);
+    MD5STEP(F2, c, d, a, b, in[3] + 0xf4d50d87, 14);
+    MD5STEP(F2, b, c, d, a, in[8] + 0x455a14ed, 20);
+    MD5STEP(F2, a, b, c, d, in[13] + 0xa9e3e905, 5);
+    MD5STEP(F2, d, a, b, c, in[2] + 0xfcefa3f8, 9);
+    MD5STEP(F2, c, d, a, b, in[7] + 0x676f02d9, 14);
+    MD5STEP(F2, b, c, d, a, in[12] + 0x8d2a4c8a, 20);
+
+    MD5STEP(F3, a, b, c, d, in[5] + 0xfffa3942, 4);
+    MD5STEP(F3, d, a, b, c, in[8] + 0x8771f681, 11);
+    MD5STEP(F3, c, d, a, b, in[11] + 0x6d9d6122, 16);
+    MD5STEP(F3, b, c, d, a, in[14] + 0xfde5380c, 23);
+    MD5STEP(F3, a, b, c, d, in[1] + 0xa4beea44, 4);
+    MD5STEP(F3, d, a, b, c, in[4] + 0x4bdecfa9, 11);
+    MD5STEP(F3, c, d, a, b, in[7] + 0xf6bb4b60, 16);
+    MD5STEP(F3, b, c, d, a, in[10] + 0xbebfbc70, 23);
+    MD5STEP(F3, a, b, c, d, in[13] + 0x289b7ec6, 4);
+    MD5STEP(F3, d, a, b, c, in[0] + 0xeaa127fa, 11);
+    MD5STEP(F3, c, d, a, b, in[3] + 0xd4ef3085, 16);
+    MD5STEP(F3, b, c, d, a, in[6] + 0x04881d05, 23);
+    MD5STEP(F3, a, b, c, d, in[9] + 0xd9d4d039, 4);
+    MD5STEP(F3, d, a, b, c, in[12] + 0xe6db99e5, 11);
+    MD5STEP(F3, c, d, a, b, in[15] + 0x1fa27cf8, 16);
+    MD5STEP(F3, b, c, d, a, in[2] + 0xc4ac5665, 23);
+
+    MD5STEP(F4, a, b, c, d, in[0] + 0xf4292244, 6);
+    MD5STEP(F4, d, a, b, c, in[7] + 0x432aff97, 10);
+    MD5STEP(F4, c, d, a, b, in[14] + 0xab9423a7, 15);
+    MD5STEP(F4, b, c, d, a, in[5] + 0xfc93a039, 21);
+    MD5STEP(F4, a, b, c, d, in[12] + 0x655b59c3, 6);
+    MD5STEP(F4, d, a, b, c, in[3] + 0x8f0ccc92, 10);
+    MD5STEP(F4, c, d, a, b, in[10] + 0xffeff47d, 15);
+    MD5STEP(F4, b, c, d, a, in[1] + 0x85845dd1, 21);
+    MD5STEP(F4, a, b, c, d, in[8] + 0x6fa87e4f, 6);
+    MD5STEP(F4, d, a, b, c, in[15] + 0xfe2ce6e0, 10);
+    MD5STEP(F4, c, d, a, b, in[6] + 0xa3014314, 15);
+    MD5STEP(F4, b, c, d, a, in[13] + 0x4e0811a1, 21);
+    MD5STEP(F4, a, b, c, d, in[4] + 0xf7537e82, 6);
+    MD5STEP(F4, d, a, b, c, in[11] + 0xbd3af235, 10);
+    MD5STEP(F4, c, d, a, b, in[2] + 0x2ad7d2bb, 15);
+    MD5STEP(F4, b, c, d, a, in[9] + 0xeb86d391, 21);
+
+    buf[0] += a;
+    buf[1] += b;
+    buf[2] += c;
+    buf[3] += d;

 }
 
 /* lh3: the following code is added by me */
 }
 
 /* lh3: the following code is added by me */


@@ -278,7 +267,7 @@ UINT4 *in;
 
 static void md5_one(const char *fn)
 {
 
 static void md5_one(const char *fn)
 {

-       unsigned char buf[4096];
+       unsigned char buf[4096], digest[16];

        MD5_CTX md5;
        int l;
        FILE *fp;
        MD5_CTX md5;
        int l;
        FILE *fp;


@@ -291,10 +280,10 @@ static void md5_one(const char *fn)
        MD5Init(&md5);
        while ((l = fread(buf, 1, 4096, fp)) > 0)
                MD5Update(&md5, buf, l);
        MD5Init(&md5);
        while ((l = fread(buf, 1, 4096, fp)) > 0)
                MD5Update(&md5, buf, l);

-       MD5Final(&md5);
+       MD5Final(digest, &md5);

        if (fp != stdin) fclose(fp);
        for (l = 0; l < 16; ++l)
        if (fp != stdin) fclose(fp);
        for (l = 0; l < 16; ++l)

-               printf("%c%c", HEX_STR[md5.digest[l]>>4&0xf], HEX_STR[md5.digest[l]&0xf]);
+               printf("%c%c", HEX_STR[digest[l]>>4&0xf], HEX_STR[digest[l]&0xf]);

        printf("  %s\n", fn);
 }
 int main(int argc, char *argv[])
        printf("  %s\n", fn);
 }
 int main(int argc, char *argv[])

diff --git a/misc/md5.h b/misc/md5.h
index 678ac277e6566d0098ea4cc35ffead3c8cb08e41..44121e4b14f7ee58fa1826e47ede4b6c12b639a8 100644 (file)
--- a/misc/md5.h
+++ b/misc/md5.h
@@ -1,68 +1,57 @@
 /*
 /*

- **********************************************************************
- ** md5.h -- Header file for implementation of MD5                   **
- ** RSA Data Security, Inc. MD5 Message Digest Algorithm             **
- ** Created: 2/17/90 RLR                                             **
- ** Revised: 12/27/90 SRD,AJ,BSK,JT Reference C version              **
- ** Revised (for MD5): RLR 4/27/91                                   **
- **   -- G modified to have y&~z instead of y&z                      **
- **   -- FF, GG, HH modified to add in last register done            **
- **   -- Access pattern: round 2 works mod 5, round 3 works mod 3    **
- **   -- distinct additive constant for each step                    **
- **   -- round 4 added, working mod 7                                **
- **********************************************************************
- */
+  This file is adapted from a program in this page:

 
 

-/*
- **********************************************************************
- ** Copyright (C) 1990, RSA Data Security, Inc. All rights reserved. **
- **                                                                  **
- ** License to copy and use this software is granted provided that   **
- ** it is identified as the "RSA Data Security, Inc. MD5 Message     **
- ** Digest Algorithm" in all material mentioning or referencing this **
- ** software or this function.                                       **
- **                                                                  **
- ** License is also granted to make and use derivative works         **
- ** provided that such works are identified as "derived from the RSA **
- ** Data Security, Inc. MD5 Message Digest Algorithm" in all         **
- ** material mentioning or referencing the derived work.             **
- **                                                                  **
- ** RSA Data Security, Inc. makes no representations concerning      **
- ** either the merchantability of this software or the suitability   **
- ** of this software for any particular purpose.  It is provided "as **
- ** is" without express or implied warranty of any kind.             **
- **                                                                  **
- ** These notices must be retained in any copies of any part of this **
- ** documentation and/or software.                                   **
- **********************************************************************
+  http://www.fourmilab.ch/md5/
+
+  The original source code does not work on 64-bit machines due to the
+  wrong typedef "uint32". I also added prototypes.
+
+  -lh3

  */
 
 #ifndef MD5_H
 #define MD5_H
 
  */
 
 #ifndef MD5_H
 #define MD5_H
 

-#include <stdint.h>
+/*  The following tests optimise behaviour on little-endian
+    machines, where there is no need to reverse the byte order
+    of 32 bit words in the MD5 computation.  By default,
+    HIGHFIRST is defined, which indicates we're running on a
+    big-endian (most significant byte first) machine, on which
+    the byteReverse function in md5.c must be invoked. However,
+    byteReverse is coded in such a way that it is an identity
+    function when run on a little-endian machine, so calling it
+    on such a platform causes no harm apart from wasting time. 
+    If the platform is known to be little-endian, we speed
+    things up by undefining HIGHFIRST, which defines
+    byteReverse as a null macro.  Doing things in this manner
+    insures we work on new platforms regardless of their byte
+    order.  */
+
+#define HIGHFIRST
+
+#if __LITTLE_ENDIAN__ != 0
+#undef HIGHFIRST
+#endif

 
 

-/* typedef a 32 bit type */
-typedef uint32_t UINT4;
+#include <stdint.h>

 
 

-/* Data structure for MD5 (Message Digest) computation */
-typedef struct {
-  UINT4 i[2];                   /* number of _bits_ handled mod 2^64 */
-  UINT4 buf[4];                                    /* scratch buffer */
-  unsigned char in[64];                              /* input buffer */
-  unsigned char digest[16];     /* actual digest after MD5Final call */
-} MD5_CTX;
+struct MD5Context {
+        uint32_t buf[4];
+        uint32_t bits[2];
+        unsigned char in[64];
+};

 
 

-#ifdef __cplusplus
-extern "C" {
-#endif
+void MD5Init(struct MD5Context *ctx);
+void MD5Update(struct MD5Context *ctx, unsigned char *buf, unsigned len);
+void MD5Final(unsigned char digest[16], struct MD5Context *ctx);

 
 

-       void MD5Init(MD5_CTX *mdContext);
-       void MD5Update(MD5_CTX *mdContext, unsigned char *inBuf, unsigned intinLen);
-       void MD5Final(MD5_CTX *mdContext);
+/*
+ * This is needed to make RSAREF happy on some MS-DOS compilers.
+ */
+typedef struct MD5Context MD5_CTX;

 
 

-#ifdef __cplusplus
-}
-#endif
+/*  Define CHECK_HARDWARE_PROPERTIES to have main,c verify
+    byte order and uint32_t settings.  */
+#define CHECK_HARDWARE_PROPERTIES

 
 

-#endif
+#endif /* !MD5_H */

diff --git a/misc/md5fa.c b/misc/md5fa.c
index c41db2d78c8b7d97f9c18d36dafeb21e9e2e8245..7a165bf74c86158ddb5e1a35799a8c4b0ffc71e5 100644 (file)
--- a/misc/md5fa.c
+++ b/misc/md5fa.c
@@ -13,7 +13,7 @@ static void md5_one(const char *fn)
        int l, i, k;
        gzFile fp;
        kseq_t *seq;
        int l, i, k;
        gzFile fp;
        kseq_t *seq;

-       unsigned char unordered[16];
+       unsigned char unordered[16], digest[16];

 
        for (l = 0; l < 16; ++l) unordered[l] = 0;
        fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
 
        for (l = 0; l < 16; ++l) unordered[l] = 0;
        fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");


@@ -31,18 +31,18 @@ static void md5_one(const char *fn)
                }
                MD5Init(&md5_one);
                MD5Update(&md5_one, (unsigned char*)seq->seq.s, k);
                }
                MD5Init(&md5_one);
                MD5Update(&md5_one, (unsigned char*)seq->seq.s, k);

-               MD5Final(&md5_one);
+               MD5Final(digest, &md5_one);

                for (l = 0; l < 16; ++l) {
                for (l = 0; l < 16; ++l) {

-                       printf("%c%c", HEX_STR[md5_one.digest[l]>>4&0xf], HEX_STR[md5_one.digest[l]&0xf]);
-                       unordered[l] ^= md5_one.digest[l];
+                       printf("%c%c", HEX_STR[digest[l]>>4&0xf], HEX_STR[digest[l]&0xf]);
+                       unordered[l] ^= digest[l];

                }
                printf("  %s  %s\n", fn, seq->name.s);
                MD5Update(&md5_all, (unsigned char*)seq->seq.s, k);
        }
                }
                printf("  %s  %s\n", fn, seq->name.s);
                MD5Update(&md5_all, (unsigned char*)seq->seq.s, k);
        }

-       MD5Final(&md5_all);
+       MD5Final(digest, &md5_all);

        kseq_destroy(seq);
        for (l = 0; l < 16; ++l)
        kseq_destroy(seq);
        for (l = 0; l < 16; ++l)

-               printf("%c%c", HEX_STR[md5_all.digest[l]>>4&0xf], HEX_STR[md5_all.digest[l]&0xf]);
+               printf("%c%c", HEX_STR[digest[l]>>4&0xf], HEX_STR[digest[l]&0xf]);

        printf("  %s  >ordered\n", fn);
        for (l = 0; l < 16; ++l)
                printf("%c%c", HEX_STR[unordered[l]>>4&0xf], HEX_STR[unordered[l]&0xf]);
        printf("  %s  >ordered\n", fn);
        for (l = 0; l < 16; ++l)
                printf("%c%c", HEX_STR[unordered[l]>>4&0xf], HEX_STR[unordered[l]&0xf]);

diff --git a/misc/psl2sam.pl b/misc/psl2sam.pl
new file mode 100755 (executable)
index 0000000..a96a6de
--- /dev/null
+++ b/misc/psl2sam.pl
@@ -0,0 +1,65 @@
+#!/usr/bin/perl -w
+
+# Author: lh3
+
+# This script calculates a score using the BLAST scoring
+# system. However, I am not sure how to count gap opens and gap
+# extensions. It seems to me that column 5-8 are not what I am
+# after. This script counts gaps from the last three columns. It does
+# not generate reference skip (N) in the CIGAR as it is not easy to
+# directly tell which gaps correspond to introns.
+
+use strict;
+use warnings;
+use Getopt::Std;
+
+my %opts = (a=>1, b=>3, q=>5, r=>2);
+getopts('a:b:q:r:', \%opts);
+die("Usage: psl2sam.pl [-a $opts{a}] [-b $opts{b}] [-q $opts{q}] [-r $opts{r}] <in.psl>\n") if (@ARGV == 0 && -t STDIN);
+
+my @stack;
+my $last = '';
+my ($a, $b, $q, $r) = ($opts{a}, $opts{b}, $opts{q}, $opts{r});
+while (<>) {
+  next unless (/^\d/);
+  my @t = split;
+  my @s;
+  my $cigar = '';
+  if ($t[8] eq '-') {
+       my $tmp = $t[11];
+       $t[11] = $t[10] - $t[12];
+       $t[12] = $t[10] - $tmp;
+  }
+  @s[0..4] = ($t[9], (($t[8] eq '+')? 0 : 16), $t[13], $t[15]+1, 0);
+  @s[6..10] = ('*', 0, 0, '*', '*');
+  $cigar .= $t[11].'H' if ($t[11]); # 5'-end clipping
+  my @x = split(',', $t[18]);
+  my @y = split(',', $t[19]);
+  my @z = split(',', $t[20]);
+  my ($y0, $z0) = ($y[0], $z[0]);
+  my ($gap_open, $gap_ext) = (0, 0, 0);
+  for (1 .. $t[17]-1) {
+       my $ly = $y[$_] - $y[$_-1] - $x[$_-1];
+       my $lz = $z[$_] - $z[$_-1] - $x[$_-1];
+       if ($ly < $lz) { # del: the reference gap is longer
+         ++$gap_open;
+         $gap_ext += $lz - $ly;
+         $cigar .= ($y[$_] - $y0) . 'M';
+         $cigar .= ($lz - $ly) . 'D';
+         ($y0, $z0) = ($y[$_], $z[$_]);
+       } elsif ($lz < $ly) { # ins: the query gap is longer
+         ++$gap_open;
+         $gap_ext += $ly - $lz;
+         $cigar .= ($z[$_] - $z0) . 'M';
+         $cigar .= ($ly - $lz) . 'I';
+         ($y0, $z0) = ($y[$_], $z[$_]);
+       }
+  }
+  $cigar .= ($t[12] - $y0) . 'M';
+  $cigar .= ($t[10] - $t[12]).'H' if ($t[10] != $t[12]); # 3'-end clipping
+  $s[5] = $cigar;
+  my $score = $a * $t[0] - $b * $t[1] - $q * $gap_open - $r * $gap_ext;
+  $score = 0 if ($score < 0);
+  $s[11] = "AS:i:$score";
+  print join("\t", @s), "\n";
+}

diff --git a/misc/samtools.pl b/misc/samtools.pl
index c014c52bcf0880d39ba23f567a31396017220142..86b285c5fa6c144a1d60ed61cb48a23bde7f134f 100755 (executable)
--- a/misc/samtools.pl
+++ b/misc/samtools.pl
@@ -6,11 +6,12 @@ use strict;
 use warnings;
 use Getopt::Std;
 
 use warnings;
 use Getopt::Std;
 

-my $version = '0.3.2 (r321)';
+my $version = '0.3.3';

 &usage if (@ARGV < 1);
 
 my $command = shift(@ARGV);
 &usage if (@ARGV < 1);
 
 my $command = shift(@ARGV);

-my %func = (showALEN=>\&showALEN, pileup2fq=>\&pileup2fq, varFilter=>\&varFilter);
+my %func = (showALEN=>\&showALEN, pileup2fq=>\&pileup2fq, varFilter=>\&varFilter,
+                       unique=>\&unique, uniqcmp=>\&uniqcmp);

 
 die("Unknown command \"$command\".\n") if (!defined($func{$command}));
 &{$func{$command}};
 
 die("Unknown command \"$command\".\n") if (!defined($func{$command}));
 &{$func{$command}};


@@ -24,9 +25,10 @@ sub showALEN {
   die(qq/Usage: samtools.pl showALEN <in.sam>\n/) if (@ARGV == 0 && -t STDIN);
   while (<>) {
        my @t = split;
   die(qq/Usage: samtools.pl showALEN <in.sam>\n/) if (@ARGV == 0 && -t STDIN);
   while (<>) {
        my @t = split;

+       next if (/^\@/ || @t < 11);

        my $l = 0;
        $_ = $t[5];
        my $l = 0;
        $_ = $t[5];

-       s/(\d+)[SMI]/$l+=$1/eg;
+       s/(\d+)[MI]/$l+=$1/eg;

        print join("\t", @t[0..5]), "\t$l\t", join("\t", @t[6..$#t]), "\n";
   }
 }
        print join("\t", @t[0..5]), "\t$l\t", join("\t", @t[6..$#t]), "\n";
   }
 }


@@ -37,7 +39,7 @@ sub showALEN {
 
 sub varFilter {
   my %opts = (d=>3, D=>100, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef);
 
 sub varFilter {
   my %opts = (d=>3, D=>100, l=>30, Q=>25, q=>10, G=>25, s=>100, w=>10, W=>10, N=>2, p=>undef);

-  getopts('pd:D:l:Q:w:W:N:G:', \%opts);
+  getopts('pq:d:D:l:Q:w:W:N:G:', \%opts);

   die(qq/
 Usage:   samtools.pl varFilter [options] <in.cns-pileup>
 
   die(qq/
 Usage:   samtools.pl varFilter [options] <in.cns-pileup>
 


@@ -65,7 +67,7 @@ Options: -Q INT    minimum RMS mapping quality for SNPs [$opts{Q}]
   my @staging; # (indel_filtering_score, flt_tag)
   while (<>) {
        my @t = split;
   my @staging; # (indel_filtering_score, flt_tag)
   while (<>) {
        my @t = split;

-       next if ($t[2] eq $t[3] || $t[3] eq '*/*'); # skip non-var sites
+       next if (uc($t[2]) eq uc($t[3]) || $t[3] eq '*/*'); # skip non-var sites

        # clear the out-of-range elements
        while (@staging) {
          last if ($staging[0][2] eq $t[0] && $staging[0][3] + $max_dist >= $t[1]);
        # clear the out-of-range elements
        while (@staging) {
          last if ($staging[0][2] eq $t[0] && $staging[0][3] + $max_dist >= $t[1]);


@@ -215,27 +217,128 @@ sub p2q_print_str {
 }
 
 #
 }
 
 #

-# varStats
+# unique

 #
 
 #
 

-sub varStats {
-  my %opts = (d=>'', c=>5);
-  getopts('d:c:', \%opts);
-  die("Usage: samtools.pl varStats [-d dbSNP.snp] [-c $opts{c}] <in.plp.snp>\n") if (@ARGV == 0 && -t STDIN);
-  my (@cnt, %hash);
-  my $col = $opts{c} - 1;
+sub unique {
+  my %opts = (f=>250.0, q=>5, r=>2, a=>1, b=>3);
+  getopts('Qf:q:r:a:b:', \%opts);
+  die("Usage: samtools.pl unique [-f $opts{f}] <in.sam>\n") if (@ARGV == 0 && -t STDIN);
+  my $last = '';
+  my $recal_Q = !defined($opts{Q});
+  my @a;

   while (<>) {
   while (<>) {

-       my @t = split;
-       if ($t[2] eq '*') {
+       my $score = -1;
+       print $_ if (/^\@/);
+       $score = $1 if (/AS:i:(\d+)/);
+       my @t = split("\t");
+       next if (@t < 11);
+       if ($score < 0) { # AS tag is unavailable
+         my $cigar = $t[5];
+         my ($mm, $go, $ge) = (0, 0, 0);
+         $cigar =~ s/(\d+)[ID]/++$go,$ge+=$1/eg;
+         $cigar = $t[5];
+         $cigar =~ s/(\d+)M/$mm+=$1/eg;
+         $score = $mm * $opts{a} - $go * $opts{q} - $ge * $opts{r}; # no mismatches...
+       }
+       $score = 1 if ($score < 1);
+       if ($t[0] ne $last) {
+         &unique_aux(\@a, $opts{f}, $recal_Q) if (@a);
+         $last = $t[0];
+       }
+       push(@a, [$score, \@t]);
+  }
+  &unique_aux(\@a, $opts{f}, $recal_Q) if (@a);
+}
+
+sub unique_aux {
+  my ($a, $fac, $is_recal) = @_;
+  my ($max, $max2, $max_i) = (0, 0, -1);
+  for (my $i = 0; $i < @$a; ++$i) {
+       if ($a->[$i][0] > $max) {
+         $max2 = $max; $max = $a->[$i][0]; $max_i = $i;
+       } elsif ($a->[$i][0] > $max2) {
+         $max2 = $a->[$i][0];
+       }
+  }
+  if ($is_recal) {
+       my $q = int($fac * ($max - $max2) / $max + .499);
+       $q = 250 if ($q > 250);
+       $a->[$max_i][1][4] = $q < 250? $q : 250;
+  }
+  print join("\t", @{$a->[$max_i][1]});
+  @$a = ();
+}
+
+#
+# uniqcmp: compare two SAM files
+#
+
+sub uniqcmp {
+  my %opts = (q=>10, s=>100);
+  getopts('pq:s:', \%opts);
+  die("Usage: samtools.pl uniqcmp <in1.sam> <in2.sam>\n") if (@ARGV < 2);
+  my ($fh, %a);
+  warn("[uniqcmp] read the first file...\n");
+  &uniqcmp_aux($ARGV[0], \%a, 0);
+  warn("[uniqcmp] read the second file...\n");
+  &uniqcmp_aux($ARGV[1], \%a, 1);
+  warn("[uniqcmp] stats...\n");
+  my @cnt;
+  $cnt[$_] = 0 for (0..9);
+  for my $x (keys %a) {
+       my $p = $a{$x};
+       my $z;
+       if (defined($p->[0]) && defined($p->[1])) {
+         $z = ($p->[0][0] == $p->[1][0] && $p->[0][1] eq $p->[1][1] && abs($p->[0][2] - $p->[1][2]) < $opts{s})? 0 : 1;
+         if ($p->[0][3] >= $opts{q} && $p->[1][3] >= $opts{q}) {
+               ++$cnt[$z*3+0];
+         } elsif ($p->[0][3] >= $opts{q}) {
+               ++$cnt[$z*3+1];
+         } elsif ($p->[1][3] >= $opts{q}) {
+               ++$cnt[$z*3+2];
+         }
+         print STDERR "$x\t$p->[0][1]:$p->[0][2]\t$p->[0][3]\t$p->[0][4]\t$p->[1][1]:$p->[1][2]\t$p->[1][3]\t$p->[1][4]\t",
+               $p->[0][5]-$p->[1][5], "\n" if ($z && defined($opts{p}) && ($p->[0][3] >= $opts{q} || $p->[1][3] >= $opts{q}));
+       } elsif (defined($p->[0])) {
+         ++$cnt[$p->[0][3]>=$opts{q}? 6 : 7];
+         print STDERR "$x\t$p->[0][1]:$p->[0][2]\t$p->[0][3]\t$p->[0][4]\t*\t0\t*\t",
+               $p->[0][5], "\n" if (defined($opts{p}) && $p->[0][3] >= $opts{q});

        } else {
        } else {

-         my $q = $t[$col];
-         $q = 99 if ($q > 99);
-         $q = int($q/10);
-         my $is_het = ($t[3] =~ /^[ACGT]$/)? 0 : 1;
-         ++$cnt[$q][$is_het];
-         $hash{$t[0],$t[1]} = $q;
+         print STDERR "$x\t*\t0\t*\t$p->[1][1]:$p->[1][2]\t$p->[1][3]\t$p->[1][4]\t",
+               -$p->[1][5], "\n" if (defined($opts{p}) && $p->[1][3] >= $opts{q});
+         ++$cnt[$p->[1][3]>=$opts{q}? 8 : 9];

        }
   }
        }
   }

+  print "Consistent (high, high):   $cnt[0]\n";
+  print "Consistent (high, low ):   $cnt[1]\n";
+  print "Consistent (low , high):   $cnt[2]\n";
+  print "Inconsistent (high, high): $cnt[3]\n";
+  print "Inconsistent (high, low ): $cnt[4]\n";
+  print "Inconsistent (low , high): $cnt[5]\n";
+  print "Second missing (high):     $cnt[6]\n";
+  print "Second missing (low ):     $cnt[7]\n";
+  print "First  missing (high):     $cnt[8]\n";
+  print "First  missing (low ):     $cnt[9]\n";
+}
+
+sub uniqcmp_aux {
+  my ($fn, $a, $which) = @_;
+  my $fh;
+  $fn = "samtools view $fn |" if ($fn =~ /\.bam/);
+  open($fh, $fn) || die;
+  while (<$fh>) {
+       my @t = split;
+       next if (@t < 11);
+#      my $l = ($t[5] =~ /^(\d+)S/)? $1 : 0;
+       my $l = 0;
+       my ($x, $nm) = (0, 0);
+       $nm = $1 if (/NM:i:(\d+)/);
+       $_ = $t[5];
+       s/(\d+)[MI]/$x+=$1/eg;
+       @{$a->{$t[0]}[$which]} = (($t[1]&0x10)? 1 : 0, $t[2], $t[3]-$l, $t[4], "$x:$nm", $x - 4 * $nm);
+  }
+  close($fh);

 }
 
 #
 }
 
 #

diff --git a/misc/wgsim_eval.pl b/misc/wgsim_eval.pl
index 99e2ac93b1c73fd37525bc963fcc4959a571cf9f..01038f1fab67c6585e2788eb19e2ae4a6e015ab2 100755 (executable)
--- a/misc/wgsim_eval.pl
+++ b/misc/wgsim_eval.pl
@@ -1,7 +1,7 @@
 #!/usr/bin/perl -w
 
 # Contact: lh3
 #!/usr/bin/perl -w
 
 # Contact: lh3

-# Version: 0.1.3
+# Version: 0.1.5

 
 use strict;
 use warnings;
 
 use strict;
 use warnings;


@@ -11,16 +11,18 @@ use Getopt::Std;
 exit;
 
 sub wgsim_eval {
 exit;
 
 sub wgsim_eval {

-  my %opts;
-  getopts('pc', \%opts);
-  die("Usage: wgsim_eval.pl [-pc] <in.sam>\n") if (@ARGV == 0 && -t STDIN);
+  my %opts = (g=>5);
+  getopts('pcg:', \%opts);
+  die("Usage: wgsim_eval.pl [-pc] [-g $opts{g}] <in.sam>\n") if (@ARGV == 0 && -t STDIN);

   my (@c0, @c1);
   my ($max_q, $flag) = (0, 0);
   my (@c0, @c1);
   my ($max_q, $flag) = (0, 0);

-  my $gap = 5;
+  my $gap = $opts{g};

   $flag |= 1 if (defined $opts{p});
   $flag |= 2 if (defined $opts{c});
   while (<>) {
   $flag |= 1 if (defined $opts{p});
   $flag |= 2 if (defined $opts{c});
   while (<>) {

-       my @t = split;
+       next if (/^\@/);
+       my @t = split("\t");
+       next if (@t < 11);

        my $line = $_;
        my ($q, $is_correct, $chr, $left, $rght) = (int($t[4]/10), 1, $t[2], $t[3], $t[3]);
        $max_q = $q if ($q > $max_q);
        my $line = $_;
        my ($q, $is_correct, $chr, $left, $rght) = (int($t[4]/10), 1, $t[2], $t[3], $t[3]);
        $max_q = $q if ($q > $max_q);


@@ -28,8 +30,11 @@ sub wgsim_eval {
        $_ = $t[5]; s/(\d+)[MDN]/$rght+=$1,'x'/eg;
        --$rght;
        # correct for soft clipping
        $_ = $t[5]; s/(\d+)[MDN]/$rght+=$1,'x'/eg;
        --$rght;
        # correct for soft clipping

-       $left -= $1 if (/^(\d+)S/);
-       $rght += $1 if (/(\d+)S$/);
+       my ($left0, $rght0) = ($left, $rght);
+       $left -= $1 if (/^(\d+)[SH]/);
+       $rght += $1 if (/(\d+)[SH]$/);
+       $left0 -= $1 if (/(\d+)[SH]$/);
+       $rght0 += $1 if (/^(\d+)[SH]/);

        # skip unmapped reads
        next if (($t[1]&0x4) || $chr eq '*');
        # parse read name and check
        # skip unmapped reads
        next if (($t[1]&0x4) || $chr eq '*');
        # parse read name and check


@@ -39,19 +44,19 @@ sub wgsim_eval {
          } else {
                if ($flag & 2) {
                  if (($t[1]&0x40) && !($t[1]&0x10)) { # F3, forward
          } else {
                if ($flag & 2) {
                  if (($t[1]&0x40) && !($t[1]&0x10)) { # F3, forward

-                       $is_correct = 0 if (abs($2 - $left) > $gap);
+                       $is_correct = 0 if (abs($2 - $left) > $gap && abs($2 - $left0) > $gap);

                  } elsif (($t[1]&0x40) && ($t[1]&0x10)) { # F3, reverse
                  } elsif (($t[1]&0x40) && ($t[1]&0x10)) { # F3, reverse

-                       $is_correct = 0 if (abs($3 - $rght) > $gap);
+                       $is_correct = 0 if (abs($3 - $rght) > $gap && abs($3 - $rght0) > $gap);

                  } elsif (($t[1]&0x80) && !($t[1]&0x10)) { # R3, forward
                  } elsif (($t[1]&0x80) && !($t[1]&0x10)) { # R3, forward

-                       $is_correct = 0 if (abs($3 - $left) > $gap);
+                       $is_correct = 0 if (abs($3 - $left) > $gap && abs($3 - $left0) > $gap);

                  } else { # R3, reverse
                  } else { # R3, reverse

-                       $is_correct = 0 if (abs($2 - $rght) > $gap);
+                       $is_correct = 0 if (abs($2 - $rght) > $gap && abs($3 - $rght0) > $gap);

                  }
                } else {
                  if ($t[1] & 0x10) { # reverse
                  }
                } else {
                  if ($t[1] & 0x10) { # reverse

-                       $is_correct = 0 if (abs($3 - $rght) > $gap); # in case of indels that are close to the end of a reads
+                       $is_correct = 0 if (abs($3 - $rght) > $gap && abs($3 - $rght0) > $gap); # in case of indels that are close to the end of a reads

                  } else {
                  } else {

-                       $is_correct = 0 if (abs($2 - $left) > $gap);
+                       $is_correct = 0 if (abs($2 - $left) > $gap && abs($2 - $left0) > $gap);

                  }
                }
          }
                  }
                }
          }

diff --git a/razf.c b/razf.c
index b56065b9e22f54fa4f8a9c31260970ff2a8808c7..a5e8f5161a8d0780a820556e0e6640fd9762bb12 100644 (file)
--- a/razf.c
+++ b/razf.c
@@ -136,6 +136,9 @@ static RAZF* razf_open_w(int fd)
 #else
 static RAZF* razf_open_w(int fd){
        RAZF *rz;
 #else
 static RAZF* razf_open_w(int fd){
        RAZF *rz;

+#ifdef _WIN32
+       setmode(fd, O_BINARY);
+#endif

        rz = calloc(1, sizeof(RAZF));
        rz->mode = 'w';
        rz->filedes = fd;
        rz = calloc(1, sizeof(RAZF));
        rz->mode = 'w';
        rz->filedes = fd;


@@ -311,6 +314,9 @@ static RAZF* razf_open_r(int fd, int _load_index){
        int n, is_be, ret;
        int64_t end;
        unsigned char c[] = "RAZF";
        int n, is_be, ret;
        int64_t end;
        unsigned char c[] = "RAZF";

+#ifdef _WIN32
+       setmode(fd, O_BINARY);
+#endif

        rz = calloc(1, sizeof(RAZF));
        rz->mode = 'r';
        rz->filedes = fd;
        rz = calloc(1, sizeof(RAZF));
        rz->mode = 'r';
        rz->filedes = fd;


@@ -382,26 +388,34 @@ static RAZF* razf_open_r(int fd, int _load_index){
 }
 
 RAZF* razf_dopen(int fd, const char *mode){
 }
 
 RAZF* razf_dopen(int fd, const char *mode){

-       if(strcasecmp(mode, "r") == 0) return razf_open_r(fd, 1);
-       else if(strcasecmp(mode, "w") == 0) return razf_open_w(fd);
+       if(strstr(mode, "r")) return razf_open_r(fd, 1);
+       else if(strstr(mode, "w")) return razf_open_w(fd);

        else return NULL;
 }
 
 RAZF* razf_dopen2(int fd, const char *mode)
 {
        else return NULL;
 }
 
 RAZF* razf_dopen2(int fd, const char *mode)
 {

-       if(strcasecmp(mode, "r") == 0) return razf_open_r(fd, 0);
-       else if(strcasecmp(mode, "w") == 0) return razf_open_w(fd);
+       if(strstr(mode, "r")) return razf_open_r(fd, 0);
+       else if(strstr(mode, "w")) return razf_open_w(fd);

        else return NULL;
 }
 
 static inline RAZF* _razf_open(const char *filename, const char *mode, int _load_index){
        int fd;
        RAZF *rz;
        else return NULL;
 }
 
 static inline RAZF* _razf_open(const char *filename, const char *mode, int _load_index){
        int fd;
        RAZF *rz;

-       if(strcasecmp(mode, "r") == 0){
+       if(strstr(mode, "r")){
+#ifdef _WIN32
+               fd = open(filename, O_RDONLY | O_BINARY);
+#else

                fd = open(filename, O_RDONLY);
                fd = open(filename, O_RDONLY);

+#endif

                rz = razf_open_r(fd, _load_index);
                rz = razf_open_r(fd, _load_index);

-       } else if(strcasecmp(mode, "w") == 0){
+       } else if(strstr(mode, "w")){
+#ifdef _WIN32
+               fd = open(filename, O_WRONLY | O_CREAT | O_TRUNC | O_BINARY, 0644);
+#else

                fd = open(filename, O_WRONLY | O_CREAT | O_TRUNC, 0644);
                fd = open(filename, O_WRONLY | O_CREAT | O_TRUNC, 0644);

+#endif

                rz = razf_open_w(fd);
        } else return NULL;
        return rz;
                rz = razf_open_w(fd);
        } else return NULL;
        return rz;

diff --git a/sam.c b/sam.c
index 45cb05cf3dd6f8a04acfd4dd92c8d250835855a2..a74c557bad9aa83d274ea44bdfe1a16dfbc4d5a6 100644 (file)
--- a/sam.c
+++ b/sam.c
@@ -1,4 +1,6 @@
 #include <string.h>
 #include <string.h>

+#include <unistd.h>
+#include "faidx.h"

 #include "sam.h"
 
 #define TYPE_BAM  1
 #include "sam.h"
 
 #define TYPE_BAM  1


@@ -75,6 +77,9 @@ samfile_t *samopen(const char *fn, const char *mode, const void *aux)
                        // open file
                        fp->x.tamw = strcmp(fn, "-")? fopen(fn, "w") : stdout;
                        if (fp->x.tamr == 0) goto open_err_ret;
                        // open file
                        fp->x.tamw = strcmp(fn, "-")? fopen(fn, "w") : stdout;
                        if (fp->x.tamr == 0) goto open_err_ret;

+                       if (strstr(mode, "X")) fp->type |= BAM_OFSTR<<2;
+                       else if (strstr(mode, "x")) fp->type |= BAM_OFHEX<<2;
+                       else fp->type |= BAM_OFDEC<<2;

                        // write header
                        if (strstr(mode, "h")) {
                                int i;
                        // write header
                        if (strstr(mode, "h")) {
                                int i;


@@ -126,7 +131,7 @@ int samwrite(samfile_t *fp, const bam1_t *b)
        if (fp == 0 || (fp->type & TYPE_READ)) return -1; // not open for writing
        if (fp->type & TYPE_BAM) return bam_write1(fp->x.bam, b);
        else {
        if (fp == 0 || (fp->type & TYPE_READ)) return -1; // not open for writing
        if (fp->type & TYPE_BAM) return bam_write1(fp->x.bam, b);
        else {

-               char *s = bam_format1(fp->header, b);
+               char *s = bam_format1_core(fp->header, b, fp->type>>2&3);

                int l = strlen(s);
                fputs(s, fp->x.tamw); fputc('\n', fp->x.tamw);
                free(s);
                int l = strlen(s);
                fputs(s, fp->x.tamw); fputc('\n', fp->x.tamw);
                free(s);


@@ -149,3 +154,23 @@ int sampileup(samfile_t *fp, int mask, bam_pileup_f func, void *func_data)
        bam_destroy1(b);
        return 0;
 }
        bam_destroy1(b);
        return 0;
 }

+
+char *samfaipath(const char *fn_ref)
+{
+       char *fn_list = 0;
+       if (fn_ref == 0) return 0;
+       fn_list = calloc(strlen(fn_ref) + 5, 1);
+       strcat(strcpy(fn_list, fn_ref), ".fai");
+       if (access(fn_list, R_OK) == -1) { // fn_list is unreadable
+               if (access(fn_ref, R_OK) == -1) {
+                       fprintf(stderr, "[samfaipath] fail to read file %s.\n", fn_ref);
+               } else {
+                       fprintf(stderr, "[samfaipath] build FASTA index...\n");
+                       if (fai_build(fn_ref) == -1) {
+                               fprintf(stderr, "[samfaipath] fail to build FASTA index.\n");
+                               free(fn_list); fn_list = 0;
+                       }
+               }
+       }
+       return fn_list;
+}

diff --git a/sam.h b/sam.h
index 970cf2daee7ca49a95280bcf1431ae35ba8f223f..0b87194e05670ba287c9debe6e22c4d79a54e353 100644 (file)
--- a/sam.h
+++ b/sam.h
@@ -15,7 +15,7 @@
 
 /*! @typedef
   @abstract SAM/BAM file handler
 
 /*! @typedef
   @abstract SAM/BAM file handler

-  @field  type    type of the handler; bit 1 for BAM and bit 2 for reading
+  @field  type    type of the handler; bit 1 for BAM, 2 for reading and bit 3-4 for flag format

   @field  bam   BAM file handler; valid if (type&1) == 1
   @field  tamr  SAM file handler for reading; valid if type == 2
   @field  tamw  SAM file handler for writing; valid if type == 0
   @field  bam   BAM file handler; valid if (type&1) == 1
   @field  tamr  SAM file handler for reading; valid if type == 2
   @field  tamw  SAM file handler for writing; valid if type == 0


@@ -41,16 +41,18 @@ extern "C" {
          @param fn SAM/BAM file name; "-" is recognized as stdin (for
          reading) or stdout (for writing).
 
          @param fn SAM/BAM file name; "-" is recognized as stdin (for
          reading) or stdout (for writing).
 

-         @param mode open mode /[rw](b?)(u?)(h?)/: 'r' for reading, 'w' for
-         writing, 'b' for BAM I/O, 'u' for uncompressed BAM output and 'h'
-         for outputing header in SAM. If 'b' present, it must immediately
-         follow 'r' or 'w'. Valid modes are "r", "w", "wh", "rb", "wb" and
-         "wbu" exclusively.
+         @param mode open mode /[rw](b?)(u?)(h?)([xX]?)/: 'r' for reading,
+         'w' for writing, 'b' for BAM I/O, 'u' for uncompressed BAM output,
+         'h' for outputing header in SAM, 'x' for HEX flag and 'X' for
+         string flag. If 'b' present, it must immediately follow 'r' or
+         'w'. Valid modes are "r", "w", "wh", "wx", "whx", "wX", "whX",
+         "rb", "wb" and "wbu" exclusively.

 
          @param aux auxiliary data; if mode[0]=='w', aux points to
 
          @param aux auxiliary data; if mode[0]=='w', aux points to

-         bam_header_t; if strcmp(mode, "rb")==0 and @SQ header lines in SAM
+         bam_header_t; if strcmp(mode, "rb")!=0 and @SQ header lines in SAM

          are absent, aux points the file name of the list of the reference;
          are absent, aux points the file name of the list of the reference;

-         aux is not used otherwise.
+         aux is not used otherwise. If @SQ header lines are present in SAM,
+         aux is not used, either.

 
          @return       SAM/BAM file handler
         */
 
          @return       SAM/BAM file handler
         */


@@ -87,6 +89,8 @@ extern "C" {
         */
        int sampileup(samfile_t *fp, int mask, bam_pileup_f func, void *data);
 
         */
        int sampileup(samfile_t *fp, int mask, bam_pileup_f func, void *data);
 

+       char *samfaipath(const char *fn_ref);
+

 #ifdef __cplusplus
 }
 #endif
 #ifdef __cplusplus
 }
 #endif

diff --git a/sam_view.c b/sam_view.c
index 02aee3c034cbc668b51295039d92a6adb1e74d51..113c6c437479237791cee16164d6fbc06cd7af8b 100644 (file)
--- a/sam_view.c
+++ b/sam_view.c
@@ -3,6 +3,7 @@
 #include <stdio.h>
 #include <unistd.h>
 #include "sam.h"
 #include <stdio.h>
 #include <unistd.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_min_mapQ = 0, g_flag_on = 0, g_flag_off = 0;
 static char *g_library, *g_rg;


@@ -31,17 +32,18 @@ static int view_func(const bam1_t *b, void *data)
        return 0;
 }
 
        return 0;
 }
 

-static int usage(void);
+static int usage(int is_long_help);

 
 int main_samview(int argc, char *argv[])
 {
        int c, is_header = 0, is_header_only = 0, is_bamin = 1, ret = 0, is_uncompressed = 0, is_bamout = 0;
 
 int 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 of_type = BAM_OFDEC, is_long_help = 0;

        samfile_t *in = 0, *out = 0;
        samfile_t *in = 0, *out = 0;

-       char in_mode[5], out_mode[5], *fn_out = 0, *fn_list = 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");
 
        /* parse command-line options */
        strcpy(in_mode, "r"); strcpy(out_mode, "w");

-       while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:ul:r:")) >= 0) {
+       while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:ul:r:xX?T:")) >= 0) {

                switch (c) {
                case 'S': is_bamin = 0; break;
                case 'b': is_bamout = 1; break;
                switch (c) {
                case 'S': is_bamin = 0; break;
                case 'b': is_bamout = 1; break;


@@ -55,22 +57,36 @@ int main_samview(int argc, char *argv[])
                case 'u': is_uncompressed = 1; break;
                case 'l': g_library = strdup(optarg); break;
                case 'r': g_rg = strdup(optarg); break;
                case 'u': is_uncompressed = 1; break;
                case 'l': g_library = strdup(optarg); break;
                case 'r': g_rg = strdup(optarg); break;

-               default: return usage();
+               case 'x': of_type = BAM_OFHEX; break;
+               case 'X': of_type = BAM_OFSTR; break;
+               case '?': is_long_help = 1; break;
+               case 'T': fn_ref = strdup(optarg); is_bamin = 0; break;
+               default: return usage(is_long_help);

                }
        }
        if (is_uncompressed) is_bamout = 1;
        if (is_header_only) is_header = 1;
        if (is_bamout) strcat(out_mode, "b");
                }
        }
        if (is_uncompressed) is_bamout = 1;
        if (is_header_only) is_header = 1;
        if (is_bamout) strcat(out_mode, "b");

+       else {
+               if (of_type == BAM_OFHEX) strcat(out_mode, "x");
+               else if (of_type == BAM_OFSTR) strcat(out_mode, "X");
+       }

        if (is_bamin) strcat(in_mode, "b");
        if (is_header) strcat(out_mode, "h");
        if (is_uncompressed) strcat(out_mode, "u");
        if (is_bamin) strcat(in_mode, "b");
        if (is_header) strcat(out_mode, "h");
        if (is_uncompressed) strcat(out_mode, "u");

-       if (argc == optind) return usage();
+       if (argc == optind) return usage(is_long_help);

 
 

+       // generate the fn_list if necessary
+       if (fn_list == 0 && fn_ref) fn_list = samfaipath(fn_ref);

        // open file handlers
        if ((in = samopen(argv[optind], in_mode, fn_list)) == 0) {
                fprintf(stderr, "[main_samview] fail to open file for reading.\n");
                goto view_end;
        }
        // open file handlers
        if ((in = samopen(argv[optind], in_mode, fn_list)) == 0) {
                fprintf(stderr, "[main_samview] fail to open file for reading.\n");
                goto view_end;
        }

+       if (in->header == 0) {
+               fprintf(stderr, "[main_samview] fail to read the header.\n");
+               goto view_end;
+       }

        if ((out = samopen(fn_out? fn_out : "-", out_mode, in->header)) == 0) {
                fprintf(stderr, "[main_samview] fail to open file for writing.\n");
                goto view_end;
        if ((out = samopen(fn_out? fn_out : "-", out_mode, in->header)) == 0) {
                fprintf(stderr, "[main_samview] fail to open file for writing.\n");
                goto view_end;


@@ -108,13 +124,13 @@ int main_samview(int argc, char *argv[])
 
 view_end:
        // close files, free and return
 
 view_end:
        // close files, free and return

-       free(fn_list); free(fn_out); free(g_library); free(g_rg);
+       free(fn_list); free(fn_ref); free(fn_out); free(g_library); free(g_rg);

        samclose(in);
        samclose(out);
        return ret;
 }
 
        samclose(in);
        samclose(out);
        return ret;
 }
 

-static int usage()
+static int usage(int is_long_help)

 {
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage:   samtools view [options] <in.bam>|<in.sam> [region1 [...]]\n\n");
 {
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage:   samtools view [options] <in.bam>|<in.sam> [region1 [...]]\n\n");


@@ -123,15 +139,20 @@ static int usage()
        fprintf(stderr, "         -H       print header only (no alignments)\n");
        fprintf(stderr, "         -S       input is SAM\n");
        fprintf(stderr, "         -u       uncompressed BAM output (force -b)\n");
        fprintf(stderr, "         -H       print header only (no alignments)\n");
        fprintf(stderr, "         -S       input is SAM\n");
        fprintf(stderr, "         -u       uncompressed BAM output (force -b)\n");

+       fprintf(stderr, "         -x       output FLAG in HEX (samtools-C specific)\n");
+       fprintf(stderr, "         -X       output FLAG in string (samtools-C specific)\n");

        fprintf(stderr, "         -t FILE  list of reference names and lengths (force -S) [null]\n");
        fprintf(stderr, "         -t FILE  list of reference names and lengths (force -S) [null]\n");

+       fprintf(stderr, "         -T FILE  reference sequence file (force -S) [null]\n");

        fprintf(stderr, "         -o FILE  output file name [stdout]\n");
        fprintf(stderr, "         -f INT   required flag, 0 for unset [0]\n");
        fprintf(stderr, "         -F INT   filtering flag, 0 for unset [0]\n");
        fprintf(stderr, "         -q INT   minimum mapping quality [0]\n");
        fprintf(stderr, "         -l STR   only output reads in library STR [null]\n");
        fprintf(stderr, "         -r STR   only output reads in read group STR [null]\n");
        fprintf(stderr, "         -o FILE  output file name [stdout]\n");
        fprintf(stderr, "         -f INT   required flag, 0 for unset [0]\n");
        fprintf(stderr, "         -F INT   filtering flag, 0 for unset [0]\n");
        fprintf(stderr, "         -q INT   minimum mapping quality [0]\n");
        fprintf(stderr, "         -l STR   only output reads in library STR [null]\n");
        fprintf(stderr, "         -r STR   only output reads in read group STR [null]\n");

-       fprintf(stderr, "\n\
-Notes:\n\
+       fprintf(stderr, "         -?       longer help\n");
+       fprintf(stderr, "\n");
+       if (is_long_help)
+               fprintf(stderr, "Notes:\n\

 \n\
   1. By default, this command assumes the file on the command line is in\n\
      the BAM format and it prints the alignments in SAM. If `-t' is\n\
 \n\
   1. By default, this command assumes the file on the command line is in\n\
      the BAM format and it prints the alignments in SAM. If `-t' is\n\


@@ -141,7 +162,7 @@ Notes:\n\
      corresponding sequence length. The `.fai' file generated by `faidx'\n\
      can be used here. This file may be empty if reads are unaligned.\n\
 \n\
      corresponding sequence length. The `.fai' file generated by `faidx'\n\
      can be used here. This file may be empty if reads are unaligned.\n\
 \n\

-  2. SAM->BAM conversion: `samtools view -bt ref.fa.fai in.sam.gz'.\n\
+  2. SAM->BAM conversion: `samtools view -bT ref.fa in.sam.gz'.\n\

 \n\
   3. BAM->SAM conversion: `samtools view in.bam'.\n\
 \n\
 \n\
   3. BAM->SAM conversion: `samtools view in.bam'.\n\
 \n\


@@ -151,6 +172,14 @@ Notes:\n\
 \n\
   5. Option `-u' is preferred over `-b' when the output is piped to\n\
      another samtools command.\n\
 \n\
   5. Option `-u' is preferred over `-b' when the output is piped to\n\
      another samtools command.\n\

+\n\
+  6. In a string FLAG, each character represents one bit with\n\
+     p=0x1 (paired), P=0x2 (properly paired), u=0x4 (unmapped),\n\
+     U=0x8 (mate unmapped), r=0x10 (reverse), R=0x20 (mate reverse)\n\
+     1=0x40 (first), 2=0x80 (second), s=0x100 (not primary), \n\
+     f=0x200 (failure) and d=0x400 (duplicate). Note that `-x' and\n\
+     `-X' are samtools-C specific. Picard and older samtools do not\n\
+     support HEX or string flags.\n\

 \n");
        return 1;
 }
 \n");
        return 1;
 }

diff --git a/samtools.1 b/samtools.1
index 45e16123948c1e72ba7193a68500e4430cc0d1c4..d2c78f1b32aeb00f4487f7944122ef9c7f5d596d 100644 (file)
--- a/samtools.1
+++ b/samtools.1
@@ -1,4 +1,4 @@
-.TH samtools 1 "6 July 2009" "samtools-0.1.5" "Bioinformatics tools"
+.TH samtools 1 "2 September 2009" "samtools-0.1.6" "Bioinformatics tools"

 .SH NAME
 .PP
 samtools - Utilities for the Sequence Alignment/Map (SAM) format
 .SH NAME
 .PP
 samtools - Utilities for the Sequence Alignment/Map (SAM) format


@@ -33,12 +33,11 @@ output (stdout). Several commands can thus be combined with Unix
 pipes. Samtools always output warning and error messages to the standard
 error output (stderr).
 
 pipes. Samtools always output warning and error messages to the standard
 error output (stderr).
 

-Samtools is also able to open a BAM (not SAM) file on a remote FTP
-server if the BAM file name starts with `ftp://'.  Samtools checks the
-current working directory for the index file and will download the index
-upon absence. Samtools achieves random FTP file access with the `REST'
-ftp command. It does not retrieve the entire alignment file unless it is
-asked to do so.
+Samtools is also able to open a BAM (not SAM) file on a remote FTP or
+HTTP server if the BAM file name starts with `ftp://' or `http://'.
+Samtools checks the current working directory for the index file and
+will download the index upon absence. Samtools does not retrieve the
+entire alignment file unless it is asked to do so.

 
 .SH COMMANDS AND OPTIONS
 
 
 .SH COMMANDS AND OPTIONS
 


@@ -73,17 +72,34 @@ Approximately the maximum required memory. [500000000]
 
 .TP
 .B merge
 
 .TP
 .B merge

-samtools merge [-n] <out.bam> <in1.bam> <in2.bam> [...]
-
-Merge multiple sorted alignments. The header of
-.I <in1.bam>
+samtools merge [-h inh.sam] [-n] <out.bam> <in1.bam> <in2.bam> [...]
+
+Merge multiple sorted alignments.
+The header reference lists of all the input BAM files, and the @SQ headers of
+.IR inh.sam ,
+if any, must all refer to the same set of reference sequences.
+The header reference list and (unless overridden by
+.BR -h )
+`@' headers of
+.I in1.bam

 will be copied to
 will be copied to

-.I <out.bam>
+.IR out.bam ,

 and the headers of other files will be ignored.
 
 .B OPTIONS:
 .RS
 .TP 8
 and the headers of other files will be ignored.
 
 .B OPTIONS:
 .RS
 .TP 8

+.B -h FILE
+Use the lines of
+.I FILE
+as `@' headers to be copied to
+.IR out.bam ,
+replacing any header lines that would otherwise be copied from
+.IR in1.bam .
+.RI ( FILE
+is actually in SAM format, though any alignment records it may contain
+are ignored.)
+.TP

 .B -n
 The input alignments are sorted by read names rather than by chromosomal
 coordinates
 .B -n
 The input alignments are sorted by read names rather than by chromosomal
 coordinates


@@ -304,9 +320,7 @@ samtools tview <in.sorted.bam> [ref.fasta]
 
 Text alignment viewer (based on the ncurses library). In the viewer,
 press `?' for help and press `g' to check the alignment start from a
 
 Text alignment viewer (based on the ncurses library). In the viewer,
 press `?' for help and press `g' to check the alignment start from a

-region in the format like `chr10:10,000,000'. Note that if the region
-showed on the screen contains no mapped reads, a blank screen will be
-seen. This is a known issue and will be improved later.
+region in the format like `chr10:10,000,000'.

 
 .RE
 
 
 .RE
 


@@ -408,6 +422,15 @@ _
 Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
 .IP o 2
 CIGAR operation P is not properly handled at the moment.
 Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
 .IP o 2
 CIGAR operation P is not properly handled at the moment.

+.IP o 2
+In merging, the input files are required to have the same number of
+reference sequences. The requirement can be relaxed. In addition,
+merging does not reconstruct the header dictionaries
+automatically. Endusers have to provide the correct header. Picard is
+better at merging.
+.IP o 2
+Samtools' rmdup does not work for single-end data and does not remove
+duplicates across chromosomes. Picard is better.

 
 .SH AUTHOR
 .PP
 
 .SH AUTHOR
 .PP


@@ -419,4 +442,4 @@ specification.
 
 .SH SEE ALSO
 .PP
 
 .SH SEE ALSO
 .PP

-Samtools website: http://samtools.sourceforge.net
+Samtools website: <http://samtools.sourceforge.net>

diff --git a/samtools.txt b/samtools.txt
new file mode 100644 (file)
index 0000000..63e6a25
--- /dev/null
+++ b/samtools.txt
@@ -0,0 +1,373 @@
+samtools(1)                  Bioinformatics tools                  samtools(1)
+
+
+
+NAME
+       samtools - Utilities for the Sequence Alignment/Map (SAM) format
+
+SYNOPSIS
+       samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
+
+       samtools sort aln.bam aln.sorted
+
+       samtools index aln.sorted.bam
+
+       samtools view aln.sorted.bam chr2:20,100,000-20,200,000
+
+       samtools merge out.bam in1.bam in2.bam in3.bam
+
+       samtools faidx ref.fasta
+
+       samtools pileup -f ref.fasta aln.sorted.bam
+
+       samtools tview aln.sorted.bam ref.fasta
+
+
+DESCRIPTION
+       Samtools  is  a  set of utilities that manipulate alignments in the BAM
+       format. It imports from and exports to the SAM (Sequence Alignment/Map)
+       format,  does  sorting,  merging  and  indexing, and allows to retrieve
+       reads in any regions swiftly.
+
+       Samtools is designed to work on a stream. It regards an input file  `-'
+       as  the  standard  input (stdin) and an output file `-' as the standard
+       output (stdout). Several commands can thus be combined with Unix pipes.
+       Samtools always output warning and error messages to the standard error
+       output (stderr).
+
+       Samtools is also able to open a BAM (not SAM) file on a remote  FTP  or
+       HTTP  server  if  the  BAM file name starts with `ftp://' or `http://'.
+       Samtools checks the current working directory for the  index  file  and
+       will  download  the  index upon absence. Samtools does not retrieve the
+       entire alignment file unless it is asked to do so.
+
+
+COMMANDS AND OPTIONS
+       import    samtools import <in.ref_list> <in.sam> <out.bam>
+
+                 Since 0.1.4, this command is an alias of:
+
+                 samtools view -bt <in.ref_list> -o <out.bam> <in.sam>
+
+
+       sort      samtools sort [-n] [-m maxMem] <in.bam> <out.prefix>
+
+                 Sort  alignments  by  leftmost  coordinates.  File  <out.pre-
+                 fix>.bam will be created. This command may also create tempo-
+                 rary files <out.prefix>.%d.bam when the whole alignment  can-
+                 not be fitted into memory (controlled by option -m).
+
+                 OPTIONS:
+
+                 -n      Sort by read names rather than by chromosomal coordi-
+                         nates
+
+                 -m INT  Approximately   the    maximum    required    memory.
+                         [500000000]
+
+
+       merge     samtools   merge   [-h   inh.sam]  [-n]  <out.bam>  <in1.bam>
+                 <in2.bam> [...]
+
+                 Merge multiple sorted alignments.  The header reference lists
+                 of  all  the input BAM files, and the @SQ headers of inh.sam,
+                 if  any,  must  all  refer  to  the  same  set  of  reference
+                 sequences.   The header reference list and (unless overridden
+                 by -h) `@' headers of in1.bam will be copied to out.bam,  and
+                 the headers of other files will be ignored.
+
+                 OPTIONS:
+
+                 -h FILE Use  the lines of FILE as `@' headers to be copied to
+                         out.bam, replacing any header lines that would other-
+                         wise  be  copied  from in1.bam.  (FILE is actually in
+                         SAM format, though any alignment records it may  con-
+                         tain are ignored.)
+
+                 -n      The  input alignments are sorted by read names rather
+                         than by chromosomal coordinates
+
+
+       index     samtools index <aln.bam>
+
+                 Index sorted alignment for fast  random  access.  Index  file
+                 <aln.bam>.bai will be created.
+
+
+       view      samtools  view  [-bhuHS]  [-t  in.refList]  [-o  output]  [-f
+                 reqFlag] [-F skipFlag] [-q minMapQ] [-l  library]  [-r  read-
+                 Group] <in.bam>|<in.sam> [region1 [...]]
+
+                 Extract/print  all or sub alignments in SAM or BAM format. If
+                 no region is specified, all the alignments will  be  printed;
+                 otherwise  only  alignments overlapping the specified regions
+                 will be output. An alignment may be given multiple  times  if
+                 it is overlapping several regions. A region can be presented,
+                 for example, in the following format: `chr2',  `chr2:1000000'
+                 or `chr2:1,000,000-2,000,000'. The coordinate is 1-based.
+
+                 OPTIONS:
+
+                 -b      Output in the BAM format.
+
+                 -u      Output uncompressed BAM. This option saves time spent
+                         on compression/decomprssion  and  is  thus  preferred
+                         when the output is piped to another samtools command.
+
+                 -h      Include the header in the output.
+
+                 -H      Output the header only.
+
+                 -S      Input is in SAM. If @SQ header lines are absent,  the
+                         `-t' option is required.
+
+                 -t FILE This  file  is  TAB-delimited. Each line must contain
+                         the reference name and the length of  the  reference,
+                         one  line  for  each  distinct  reference; additional
+                         fields are ignored. This file also defines the  order
+                         of  the  reference  sequences  in sorting. If you run
+                         `samtools faidx <ref.fa>', the resultant  index  file
+                         <ref.fa>.fai  can be used as this <in.ref_list> file.
+
+                 -o FILE Output file [stdout]
+
+                 -f INT  Only output alignments with all bits in  INT  present
+                         in the FLAG field. INT can be in hex in the format of
+                         /^0x[0-9A-F]+/ [0]
+
+                 -F INT  Skip alignments with bits present in INT [0]
+
+                 -q INT  Skip alignments with MAPQ smaller than INT [0]
+
+                 -l STR  Only output reads in library STR [null]
+
+                 -r STR  Only output reads in read group STR [null]
+
+
+       faidx     samtools faidx <ref.fasta> [region1 [...]]
+
+                 Index reference sequence in the FASTA format or extract  sub-
+                 sequence  from  indexed  reference  sequence. If no region is
+                 specified,   faidx   will   index   the   file   and   create
+                 <ref.fasta>.fai  on the disk. If regions are speficified, the
+                 subsequences will be retrieved and printed to stdout  in  the
+                 FASTA  format.  The  input file can be compressed in the RAZF
+                 format.
+
+
+       pileup    samtools  pileup  [-f  in.ref.fasta]  [-t  in.ref_list]   [-l
+                 in.site_list]    [-iscgS2]   [-T   theta]   [-N   nHap]   [-r
+                 pairDiffRate] <in.bam>|<in.sam>
+
+                 Print the alignment in the pileup format. In the pileup  for-
+                 mat,  each  line represents a genomic position, consisting of
+                 chromosome name, coordinate, reference base, read bases, read
+                 qualities  and  alignment  mapping  qualities. Information on
+                 match, mismatch, indel, strand, mapping quality and start and
+                 end  of  a  read  are all encoded at the read base column. At
+                 this column, a dot stands for a match to the  reference  base
+                 on  the  forward  strand,  a comma for a match on the reverse
+                 strand, `ACGTN' for a mismatch  on  the  forward  strand  and
+                 `acgtn'  for  a  mismatch  on  the  reverse strand. A pattern
+                 `\+[0-9]+[ACGTNacgtn]+'  indicates  there  is  an   insertion
+                 between  this reference position and the next reference posi-
+                 tion. The length of the insertion is given by the integer  in
+                 the  pattern, followed by the inserted sequence. Similarly, a
+                 pattern `-[0-9]+[ACGTNacgtn]+' represents a deletion from the
+                 reference.  The deleted bases will be presented as `*' in the
+                 following lines. Also at the read base column, a  symbol  `^'
+                 marks  the start of a read segment which is a contiguous sub-
+                 sequence on the read separated by `N/S/H'  CIGAR  operations.
+                 The  ASCII  of the character following `^' minus 33 gives the
+                 mapping quality. A symbol `$' marks the end of  a  read  seg-
+                 ment.
+
+                 If  option -c is applied, the consensus base, consensus qual-
+                 ity, SNP quality and RMS mapping quality of the reads  cover-
+                 ing  the  site  will be inserted between the `reference base'
+                 and the `read bases' columns. An indel occupies an additional
+                 line.  Each  indel  line consists of chromosome name, coordi-
+                 nate, a star, the genotype, consensus quality,  SNP  quality,
+                 RMS mapping quality, # covering reads, the first alllele, the
+                 second allele, # reads supporting the first allele,  #  reads
+                 supporting  the  second  allele and # reads containing indels
+                 different from the top two alleles.
+
+                 OPTIONS:
+
+
+                 -s        Print the mapping quality as the last column.  This
+                           option  makes  the output easier to parse, although
+                           this format is not space efficient.
+
+
+                 -S        The input file is in SAM.
+
+
+                 -i        Only output pileup lines containing indels.
+
+
+                 -f FILE   The reference sequence in the FASTA  format.  Index
+                           file FILE.fai will be created if absent.
+
+
+                 -M INT    Cap mapping quality at INT [60]
+
+
+                 -t FILE   List  of  reference  names ane sequence lengths, in
+                           the format described for  the  import  command.  If
+                           this  option is present, samtools assumes the input
+                           <in.alignment>  is  in  SAM  format;  otherwise  it
+                           assumes in BAM format.
+
+
+                 -l FILE   List  of sites at which pileup is output. This file
+                           is space  delimited.  The  first  two  columns  are
+                           required  to  be chromosome and 1-based coordinate.
+                           Additional columns are ignored. It  is  recommended
+                           to use option -s together with -l as in the default
+                           format we may not know the mapping quality.
+
+
+                 -c        Call the consensus  sequence  using  MAQ  consensus
+                           model. Options -T, -N, -I and -r are only effective
+                           when -c or -g is in use.
+
+
+                 -g        Generate genotype likelihood in  the  binary  GLFv3
+                           format. This option suppresses -c, -i and -s.
+
+
+                 -T FLOAT  The  theta parameter (error dependency coefficient)
+                           in the maq consensus calling model [0.85]
+
+
+                 -N INT    Number of haplotypes in the sample (>=2) [2]
+
+
+                 -r FLOAT  Expected fraction of differences between a pair  of
+                           haplotypes [0.001]
+
+
+                 -I INT    Phred  probability  of an indel in sequencing/prep.
+                           [40]
+
+
+
+       tview     samtools tview <in.sorted.bam> [ref.fasta]
+
+                 Text alignment viewer (based on the ncurses library). In  the
+                 viewer,  press `?' for help and press `g' to check the align-
+                 ment   start   from   a   region   in   the    format    like
+                 `chr10:10,000,000'.
+
+
+
+       fixmate   samtools fixmate <in.nameSrt.bam> <out.bam>
+
+                 Fill in mate coordinates, ISIZE and mate related flags from a
+                 name-sorted alignment.
+
+
+       rmdup     samtools rmdup <input.srt.bam> <out.bam>
+
+                 Remove potential PCR duplicates: if multiple read pairs  have
+                 identical  external  coordinates,  only  retain the pair with
+                 highest mapping quality.  This command  ONLY  works  with  FR
+                 orientation and requires ISIZE is correctly set.
+
+
+
+       rmdupse   samtools rmdupse <input.srt.bam> <out.bam>
+
+                 Remove potential duplicates for single-ended reads. This com-
+                 mand will treat all reads as single-ended even  if  they  are
+                 paired in fact.
+
+
+
+       fillmd    samtools fillmd [-e] <aln.bam> <ref.fasta>
+
+                 Generate  the  MD tag. If the MD tag is already present, this
+                 command will give a warning if the MD tag generated  is  dif-
+                 ferent from the existing tag.
+
+                 OPTIONS:
+
+                 -e      Convert  a  the  read base to = if it is identical to
+                         the aligned reference base.  Indel  caller  does  not
+                         support the = bases at the moment.
+
+
+
+SAM FORMAT
+       SAM  is  TAB-delimited.  Apart from the header lines, which are started
+       with the `@' symbol, each alignment line consists of:
+
+
+       +----+-------+----------------------------------------------------------+
+       |Col | Field |                       Description                        |
+       +----+-------+----------------------------------------------------------+
+       | 1  | QNAME | Query (pair) NAME                                        |
+       | 2  | FLAG  | bitwise FLAG                                             |
+       | 3  | RNAME | Reference sequence NAME                                  |
+       | 4  | POS   | 1-based leftmost POSition/coordinate of clipped sequence |
+       | 5  | MAPQ  | MAPping Quality (Phred-scaled)                           |
+       | 6  | CIAGR | extended CIGAR string                                    |
+       | 7  | MRNM  | Mate Reference sequence NaMe (`=' if same as RNAME)      |
+       | 8  | MPOS  | 1-based Mate POSistion                                   |
+       | 9  | ISIZE | Inferred insert SIZE                                     |
+       |10  | SEQ   | query SEQuence on the same strand as the reference       |
+       |11  | QUAL  | query QUALity (ASCII-33 gives the Phred base quality)    |
+       |12  | OPT   | variable OPTional fields in the format TAG:VTYPE:VALUE   |
+       +----+-------+----------------------------------------------------------+
+
+       Each bit in the FLAG field is defined as:
+
+
+             +-------+--------------------------------------------------+
+             | Flag  |                   Description                    |
+             +-------+--------------------------------------------------+
+             |0x0001 | the read is paired in sequencing                 |
+             |0x0002 | the read is mapped in a proper pair              |
+             |0x0004 | the query sequence itself is unmapped            |
+             |0x0008 | the mate is unmapped                             |
+             |0x0010 | strand of the query (1 for reverse)              |
+             |0x0020 | strand of the mate                               |
+             |0x0040 | the read is the first read in a pair             |
+             |0x0080 | the read is the second read in a pair            |
+             |0x0100 | the alignment is not primary                     |
+             |0x0200 | the read fails platform/vendor quality checks    |
+             |0x0400 | the read is either a PCR or an optical duplicate |
+             +-------+--------------------------------------------------+
+
+LIMITATIONS
+       o Unaligned  words  used  in  bam_import.c,  bam_endian.h,  bam.c   and
+         bam_aux.c.
+
+       o CIGAR operation P is not properly handled at the moment.
+
+       o In  merging,  the input files are required to have the same number of
+         reference sequences. The requirement can  be  relaxed.  In  addition,
+         merging  does  not reconstruct the header dictionaries automatically.
+         Endusers have to provide the correct  header.  Picard  is  better  at
+         merging.
+
+       o Samtools' rmdup does not work for single-end data and does not remove
+         duplicates across chromosomes. Picard is better.
+
+
+AUTHOR
+       Heng Li from the Sanger Institute wrote the C version of samtools.  Bob
+       Handsaker from the Broad Institute implemented the BGZF library and Jue
+       Ruan from Beijing Genomics Institute wrote the  RAZF  library.  Various
+       people  in the 1000Genomes Project contributed to the SAM format speci-
+       fication.
+
+
+SEE ALSO
+       Samtools website: <http://samtools.sourceforge.net>
+
+
+
+samtools-0.1.6                 2 September 2009                    samtools(1)

diff --git a/source.dot b/source.dot
deleted file mode 100644 (file)
index 1735774..0000000
--- a/source.dot
+++ /dev/null
@@ -1,19 +0,0 @@
-digraph {
-  faidx[label="faidx.c\n(faidx)"]
-  import[label="bam_import.c\n(import)"]
-  plcmd[label="bam_plcmd.c\n(pileup)"]
-  sort[label="bam_sort.c\n(sort, merge)"]
-  index[label="bam_index.c\n(index)"]
-  tview[label="bam_tview.c\n(tview)"]
-  glf[label="glf.c\n(glfview)"]
-  rmdup[label="bam_rmdup.c\n(rmdup)"]
-  fixmate[label="bam_mate.c\n(fixmate)"]
-  "bam_aux.c" -> {"bam.c", import}
-  glf -> {"bam_maqcns.c", plcmd}
-  "bgzf.c" -> {"bam.c", glf}
-  "bam.c" -> {index, "bam_pileup.c", sort, import, rmdup, fixmate}
-  "bam_pileup.c" -> {"bam_lpileup.c", plcmd}
-  {"bam_lpileup.c", index, faidx, "bam_maqcns.c"} -> tview
-  {import, faidx, "bam_maqcns.c"} -> plcmd
-  {tview, plcmd, faidx, sort, import, index, glf, rmdup, fixmate} -> "bamtk.c\n(view)"
-}
\ No newline at end of file