[med-svn] [libtabixpp] 01/03: Imported Upstream version 0.0.20151117
Sascha Steinbiss
sascha-guest at moszumanska.debian.org
Sun Dec 6 18:59:41 UTC 2015
This is an automated email from the git hooks/post-receive script.
sascha-guest pushed a commit to branch master
in repository libtabixpp.
commit 8af1d9f1ff7610b84993e85c293b1bc7b173026c
Author: Sascha Steinbiss <sascha at steinbiss.name>
Date: Sun Dec 6 12:48:03 2015 +0000
Imported Upstream version 0.0.20151117
---
.gitmodules | 3 +
ChangeLog | 593 ------------------------------
Makefile | 90 ++---
NEWS | 126 -------
README | 6 -
README.md | 5 +
TabixReader.java | 395 --------------------
bam_endian.h | 42 ---
bedidx.c | 156 --------
bgzf.c | 714 ------------------------------------
bgzf.h | 157 --------
bgzip.c | 206 -----------
example.gtf.gz | Bin 3778 -> 0 bytes
example.gtf.gz.tbi | Bin 259 -> 0 bytes
index.c | 998 --------------------------------------------------
khash.h | 486 ------------------------
knetfile.c | 632 --------------------------------
knetfile.h | 75 ----
kseq.h | 227 ------------
ksort.h | 271 --------------
kstring.c | 165 ---------
kstring.h | 68 ----
main.c | 290 ---------------
perl/MANIFEST | 8 -
perl/Makefile.PL | 8 -
perl/Tabix.pm | 76 ----
perl/Tabix.xs | 71 ----
perl/TabixIterator.pm | 41 ---
perl/t/01local.t | 28 --
perl/t/02remote.t | 28 --
perl/typemap | 3 -
python/setup.py | 55 ---
python/tabixmodule.c | 408 ---------------------
python/test.py | 91 -----
tabix.1 | 132 -------
tabix.cpp | 85 +++--
tabix.h | 145 --------
tabix.hpp | 17 +-
tabix.py | 87 -----
tabix.tex | 121 ------
40 files changed, 107 insertions(+), 7002 deletions(-)
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..9e410d0
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "htslib"]
+ path = htslib
+ url = https://github.com/samtools/htslib.git
diff --git a/ChangeLog b/ChangeLog
deleted file mode 100644
index fd335b8..0000000
--- a/ChangeLog
+++ /dev/null
@@ -1,593 +0,0 @@
-------------------------------------------------------------------------
-r942 | lh3lh3 | 2011-03-31 16:39:50 -0400 (Thu, 31 Mar 2011) | 2 lines
-Changed paths:
- M /trunk/tabix/main.c
-
-update version number
-
-------------------------------------------------------------------------
-r940 | lh3lh3 | 2011-03-31 16:38:03 -0400 (Thu, 31 Mar 2011) | 2 lines
-Changed paths:
- M /trunk/tabix/bedidx.c
- M /trunk/tabix/main.c
-
-fixed two bugs due to recent changes
-
-------------------------------------------------------------------------
-r939 | lh3lh3 | 2011-03-31 16:12:21 -0400 (Thu, 31 Mar 2011) | 2 lines
-Changed paths:
- M /trunk/tabix/bgzf.c
- M /trunk/tabix/bgzf.h
- M /trunk/tabix/main.c
-
-update to the latest bgzf.*
-
-------------------------------------------------------------------------
-r938 | lh3lh3 | 2011-03-31 16:02:21 -0400 (Thu, 31 Mar 2011) | 2 lines
-Changed paths:
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.h
-
-BED support
-
-------------------------------------------------------------------------
-r937 | lh3lh3 | 2011-03-31 15:03:49 -0400 (Thu, 31 Mar 2011) | 2 lines
-Changed paths:
- M /trunk/tabix/Makefile
- A /trunk/tabix/bedidx.c
- M /trunk/tabix/example.gtf.gz.tbi
- M /trunk/tabix/index.c
- A /trunk/tabix/kseq.h
- M /trunk/tabix/tabix.h
-
-restructure get_intv() for BED support
-
-------------------------------------------------------------------------
-r919 | petulda | 2011-02-24 10:14:14 -0500 (Thu, 24 Feb 2011) | 1 line
-Changed paths:
- M /trunk/tabix/bgzf.c
- M /trunk/tabix/bgzf.h
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
-New -r (reheader) option for efficient header replacement.
-------------------------------------------------------------------------
-r915 | lh3lh3 | 2011-02-22 09:50:57 -0500 (Tue, 22 Feb 2011) | 2 lines
-Changed paths:
- A /trunk/tabix/python
- A /trunk/tabix/python/setup.py (from /trunk/tabix/setup.py:914)
- A /trunk/tabix/python/tabixmodule.c (from /trunk/tabix/tabixmodule.c:914)
- A /trunk/tabix/python/test.py (from /trunk/tabix/test.py:914)
- D /trunk/tabix/setup.py
- D /trunk/tabix/tabixmodule.c
- D /trunk/tabix/test.py
-
-move to a new python/ directory
-
-------------------------------------------------------------------------
-r914 | lh3lh3 | 2011-02-22 09:49:35 -0500 (Tue, 22 Feb 2011) | 2 lines
-Changed paths:
- A /trunk/tabix/setup.py
- A /trunk/tabix/tabixmodule.c
- A /trunk/tabix/test.py
-
-CPython C-API by Hyeshik Chang
-
-------------------------------------------------------------------------
-r904 | petulda | 2011-01-28 08:06:27 -0500 (Fri, 28 Jan 2011) | 1 line
-Changed paths:
- M /trunk/tabix/index.c
-
-Check the number of fields on each line and exit nicely without segfault
-------------------------------------------------------------------------
-r901 | petulda | 2011-01-21 06:45:37 -0500 (Fri, 21 Jan 2011) | 1 line
-Changed paths:
- M /trunk/tabix/main.c
-
-Fix: Complain only when VCF is newer, not newer or same mtime
-------------------------------------------------------------------------
-r900 | petulda | 2011-01-21 04:23:04 -0500 (Fri, 21 Jan 2011) | 1 line
-Changed paths:
- M /trunk/tabix/main.c
-
-Prevent the common user mistake and check the timestamps of the vcf and index file
-------------------------------------------------------------------------
-r876 | lh3lh3 | 2010-12-08 12:38:45 -0500 (Wed, 08 Dec 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/ChangeLog
- M /trunk/tabix/NEWS
- M /trunk/tabix/main.c
-
-Release tabix-0.2.3
-
-------------------------------------------------------------------------
-r875 | lh3lh3 | 2010-12-08 12:28:35 -0500 (Wed, 08 Dec 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/ChangeLog
- M /trunk/tabix/index.c
-
-Fixed a minor bug in generating index
-
-------------------------------------------------------------------------
-r855 | petulda | 2010-11-25 11:50:13 -0500 (Thu, 25 Nov 2010) | 1 line
-Changed paths:
- M /trunk/tabix/main.c
-
-Disable "unknown target name or minus interval" warning.
-------------------------------------------------------------------------
-r775 | petulda | 2010-10-26 15:02:30 -0400 (Tue, 26 Oct 2010) | 1 line
-Changed paths:
- M /trunk/tabix/main.c
-
-Added -h option to print header lines
-------------------------------------------------------------------------
-r742 | jmarshall | 2010-09-27 06:47:23 -0400 (Mon, 27 Sep 2010) | 2 lines
-Changed paths:
- M /trunk/tabix
-
-Add svn:ignore properties for intermediate and generated files.
-
-------------------------------------------------------------------------
-r725 | lh3lh3 | 2010-09-15 13:01:53 -0400 (Wed, 15 Sep 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/bgzip.c
-
-patches by Peter Chines
-
-------------------------------------------------------------------------
-r714 | lh3lh3 | 2010-09-07 10:13:25 -0400 (Tue, 07 Sep 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/TabixReader.java
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
-fixed a bug in C/Java when n_off == 0
-
-------------------------------------------------------------------------
-r712 | lh3lh3 | 2010-09-03 09:21:23 -0400 (Fri, 03 Sep 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/TabixReader.java
-
-fixed a bug in parsing region strings
-
-------------------------------------------------------------------------
-r700 | petulda | 2010-08-25 10:42:37 -0400 (Wed, 25 Aug 2010) | 1 line
-Changed paths:
- M /trunk/tabix/main.c
-
-Fix: Exit with an error rather than segfault when index is not present and region is queried
-------------------------------------------------------------------------
-r696 | petulda | 2010-08-24 10:24:12 -0400 (Tue, 24 Aug 2010) | 1 line
-Changed paths:
- M /trunk/tabix/bgzf.c
- M /trunk/tabix/bgzf.h
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
-Complain about not-bgzipped files and check for noncontinuous chromosome blocks
-------------------------------------------------------------------------
-r603 | lh3lh3 | 2010-06-28 10:49:39 -0400 (Mon, 28 Jun 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/NEWS
- M /trunk/tabix/TabixReader.java
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
-Release tabix-0.2.2
-
-------------------------------------------------------------------------
-r597 | lh3lh3 | 2010-06-13 21:08:29 -0400 (Sun, 13 Jun 2010) | 3 lines
-Changed paths:
- M /trunk/tabix/index.c
-
-Change the namespace of sorting, to avoid function name collision with samtools.
-
-
-------------------------------------------------------------------------
-r582 | lh3lh3 | 2010-06-03 10:40:25 -0400 (Thu, 03 Jun 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/NEWS
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.py
-
-Release tabix-0.2.1
-
-------------------------------------------------------------------------
-r581 | lh3lh3 | 2010-05-24 14:24:24 -0400 (Mon, 24 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/tabix.py
-
-OOP interface with the help from Aaron Quinlan
-
-------------------------------------------------------------------------
-r580 | lh3lh3 | 2010-05-23 23:36:05 -0400 (Sun, 23 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/tabix.py
-
-minor change
-
-------------------------------------------------------------------------
-r579 | lh3lh3 | 2010-05-23 23:25:24 -0400 (Sun, 23 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/tabix.py
-
-For Snow Leopard compatibility
-
-------------------------------------------------------------------------
-r575 | lh3lh3 | 2010-05-12 19:31:27 -0400 (Wed, 12 May 2010) | 4 lines
-Changed paths:
- M /trunk/tabix/Makefile
- M /trunk/tabix/index.c
- M /trunk/tabix/tabix.h
- A /trunk/tabix/tabix.py
-
- * optionally generate shared library for Mac and Linux
- * added a python script that directly calls the shared library
- * added a new API for easy python access
-
-------------------------------------------------------------------------
-r574 | lh3lh3 | 2010-05-11 12:14:27 -0400 (Tue, 11 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/ChangeLog
- M /trunk/tabix/NEWS
- M /trunk/tabix/perl/Tabix.pm
- M /trunk/tabix/perl/TabixIterator.pm
- M /trunk/tabix/tabix.1
-
-Release tabix-0.2.0
-
-------------------------------------------------------------------------
-r573 | lh3lh3 | 2010-05-11 12:08:30 -0400 (Tue, 11 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/Makefile
-
-Added -fPIC
-
-------------------------------------------------------------------------
-r572 | lh3lh3 | 2010-05-11 11:59:07 -0400 (Tue, 11 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/perl/MANIFEST
-
-update
-
-------------------------------------------------------------------------
-r571 | lh3lh3 | 2010-05-11 11:56:54 -0400 (Tue, 11 May 2010) | 4 lines
-Changed paths:
- A /trunk/tabix/example.gtf.gz
- A /trunk/tabix/example.gtf.gz.tbi
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/perl/MANIFEST
- M /trunk/tabix/perl/Tabix.pm
- M /trunk/tabix/perl/Tabix.xs
- A /trunk/tabix/perl/TabixIterator.pm
- A /trunk/tabix/perl/t
- A /trunk/tabix/perl/t/01local.t
- A /trunk/tabix/perl/t/02remote.t
- M /trunk/tabix/tabix.1
- M /trunk/tabix/tabix.h
-
- * improved C/Perl APIs
- * added test for Perl
- * added an tiny example
-
-------------------------------------------------------------------------
-r570 | lh3lh3 | 2010-05-11 01:04:21 -0400 (Tue, 11 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/TabixReader.java
-
-fixed the same issue in java
-
-------------------------------------------------------------------------
-r569 | lh3lh3 | 2010-05-11 01:03:24 -0400 (Tue, 11 May 2010) | 3 lines
-Changed paths:
- M /trunk/tabix/index.c
- M /trunk/tabix/perl/Tabix.pm
- M /trunk/tabix/perl/Tabix.xs
-
- * fixed a potential issue in index.c
- * improve perl APIs
-
-------------------------------------------------------------------------
-r568 | lh3lh3 | 2010-05-10 23:46:21 -0400 (Mon, 10 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/perl/Tabix.xs
-
-return an array from get_names()
-
-------------------------------------------------------------------------
-r567 | lh3lh3 | 2010-05-10 23:38:46 -0400 (Mon, 10 May 2010) | 4 lines
-Changed paths:
- M /trunk/tabix/TabixReader.java
- M /trunk/tabix/index.c
- A /trunk/tabix/perl
- A /trunk/tabix/perl/MANIFEST
- A /trunk/tabix/perl/Makefile.PL
- A /trunk/tabix/perl/Tabix.pm
- A /trunk/tabix/perl/Tabix.xs
- A /trunk/tabix/perl/typemap
- M /trunk/tabix/tabix.h
-
- * added the initial perl binding. The interface needs to be improved.
- * added a new API for perl binding
- * fixed a potential bug in java.
-
-------------------------------------------------------------------------
-r565 | lh3lh3 | 2010-05-09 23:24:35 -0400 (Sun, 09 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/main.c
-
-Release tabix-0.1.6
-
-------------------------------------------------------------------------
-r564 | lh3lh3 | 2010-05-09 23:01:49 -0400 (Sun, 09 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/index.c
-
-fixed a typo
-
-------------------------------------------------------------------------
-r563 | lh3lh3 | 2010-05-09 22:58:26 -0400 (Sun, 09 May 2010) | 2 lines
-Changed paths:
- A /trunk/tabix/ChangeLog
- M /trunk/tabix/NEWS
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.h
-
-If nothing bad happens, this will become 0.1.6
-
-------------------------------------------------------------------------
-r562 | lh3lh3 | 2010-05-09 19:43:56 -0400 (Sun, 09 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/index.c
-
-Fixed a bug
-
-------------------------------------------------------------------------
-r560 | lh3lh3 | 2010-05-05 10:59:09 -0400 (Wed, 05 May 2010) | 3 lines
-Changed paths:
- A /trunk/tabix/NEWS
- M /trunk/tabix/TabixReader.java
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.1
- M /trunk/tabix/tabix.h
-
- * Release tabix-0.1.5 (r560)
- * Improve seeking efficiency. Index file needs to be rebuilt.
-
-------------------------------------------------------------------------
-r559 | lh3lh3 | 2010-05-04 23:11:42 -0400 (Tue, 04 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/main.c
-
-Release tabix-0.1.4 (r559)
-
-------------------------------------------------------------------------
-r558 | lh3lh3 | 2010-05-01 12:48:01 -0400 (Sat, 01 May 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/TabixReader.java
-
-implement SAM/VCF support; NOT tested yet
-
-------------------------------------------------------------------------
-r557 | lh3lh3 | 2010-05-01 00:42:34 -0400 (Sat, 01 May 2010) | 2 lines
-Changed paths:
- A /trunk/tabix/TabixReader.java
-
-The Java implementation of tabix.
-
-------------------------------------------------------------------------
-r556 | lh3lh3 | 2010-04-30 22:34:07 -0400 (Fri, 30 Apr 2010) | 4 lines
-Changed paths:
- M /trunk/tabix/index.c
- M /trunk/tabix/knetfile.c
- M /trunk/tabix/main.c
-
- * tabix-0.1.3-3 (r556)
- * fixed a small memory leak in knetfile
- * fixed a minor bug for remote downloading
-
-------------------------------------------------------------------------
-r555 | lh3lh3 | 2010-04-30 22:15:12 -0400 (Fri, 30 Apr 2010) | 4 lines
-Changed paths:
- M /trunk/tabix/Makefile
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
- * tabix-0.1.3-2 (r555)
- * do not overwrite index file by default
- * a little code cleanup
-
-------------------------------------------------------------------------
-r554 | lh3lh3 | 2010-04-30 21:44:31 -0400 (Fri, 30 Apr 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/index.c
-
-fixed a potential bug for UCSC-like coordinate
-
-------------------------------------------------------------------------
-r553 | lh3lh3 | 2010-04-28 17:43:41 -0400 (Wed, 28 Apr 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/tabix.tex
-
-minor clarification to the format spec
-
-------------------------------------------------------------------------
-r552 | lh3lh3 | 2010-04-28 16:33:07 -0400 (Wed, 28 Apr 2010) | 3 lines
-Changed paths:
- M /trunk/tabix/Makefile
- M /trunk/tabix/bgzip.c
- A /trunk/tabix/tabix.tex
-
- * added the format specification
- * fixed a typo in bgzip
-
-------------------------------------------------------------------------
-r550 | petulda | 2010-04-22 11:03:24 -0400 (Thu, 22 Apr 2010) | 1 line
-Changed paths:
- M /trunk/tabix/bgzip.c
-
-The behaviour changed slightly to mimic gzip. Detect if std descriptors are connected to the terminal.
-------------------------------------------------------------------------
-r549 | petulda | 2010-04-22 09:46:10 -0400 (Thu, 22 Apr 2010) | 1 line
-Changed paths:
- M /trunk/tabix/bgzip.c
-
-Fix in src/dst file detection and slight change of behaviour
-------------------------------------------------------------------------
-r548 | petulda | 2010-04-19 04:39:46 -0400 (Mon, 19 Apr 2010) | 1 line
-Changed paths:
- M /trunk/tabix/index.c
-
-Close file descriptor in ti_list_chromosomes
-------------------------------------------------------------------------
-r547 | petulda | 2010-04-16 09:27:11 -0400 (Fri, 16 Apr 2010) | 1 line
-Changed paths:
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.h
-
-Added the -l option for listing chromosomes
-------------------------------------------------------------------------
-r544 | lh3lh3 | 2010-03-29 10:58:48 -0400 (Mon, 29 Mar 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/main.c
-
-removed a line of debugging code
-
-------------------------------------------------------------------------
-r543 | lh3lh3 | 2010-03-19 12:29:16 -0400 (Fri, 19 Mar 2010) | 3 lines
-Changed paths:
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.1
-
- * tabix-0.1.3 (r543)
- * fixed another off-by-one bug
-
-------------------------------------------------------------------------
-r542 | lh3lh3 | 2010-03-16 22:35:52 -0400 (Tue, 16 Mar 2010) | 2 lines
-Changed paths:
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.1
-
-Release tabix-0.1.1
-
-------------------------------------------------------------------------
-r506 | lh3lh3 | 2009-11-02 23:20:12 -0500 (Mon, 02 Nov 2009) | 2 lines
-Changed paths:
- M /trunk/tabix/main.c
-
-Release tabix-0.1.0
-
-------------------------------------------------------------------------
-r505 | lh3lh3 | 2009-11-02 23:15:49 -0500 (Mon, 02 Nov 2009) | 2 lines
-Changed paths:
- A /trunk/tabix/tabix.1
-
-documentation
-
-------------------------------------------------------------------------
-r504 | lh3lh3 | 2009-11-02 11:08:18 -0500 (Mon, 02 Nov 2009) | 5 lines
-Changed paths:
- M /trunk/tabix/Makefile
- M /trunk/tabix/bgzip.c
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.h
-
- * tabix-0.0.0-5 (r504)
- * fixed a critical bug in fetching data (a typo in fact)
- * support SAM (tested on ex1.sam) and VCF (not tested)
- * improve the command-line interface
-
-------------------------------------------------------------------------
-r503 | lh3lh3 | 2009-11-02 10:04:43 -0500 (Mon, 02 Nov 2009) | 3 lines
-Changed paths:
- M /trunk/tabix/Makefile
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
- * tabix-0.0.0-4 (r503)
- * index files are bgzf compressed
-
-------------------------------------------------------------------------
-r502 | lh3lh3 | 2009-11-02 09:47:25 -0500 (Mon, 02 Nov 2009) | 4 lines
-Changed paths:
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
- M /trunk/tabix/tabix.h
-
- * tabix-0.0.0-3 (r502)
- * support meta lines (not tested)
- * I am going to make the index file in the BGZF format
-
-------------------------------------------------------------------------
-r501 | lh3lh3 | 2009-11-01 22:03:07 -0500 (Sun, 01 Nov 2009) | 3 lines
-Changed paths:
- M /trunk/tabix/Makefile
- M /trunk/tabix/bgzf.h
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
- * tabix-0.0.0-2 (r501)
- * accelerate ti_readline()
-
-------------------------------------------------------------------------
-r500 | lh3lh3 | 2009-11-01 20:49:52 -0500 (Sun, 01 Nov 2009) | 3 lines
-Changed paths:
- M /trunk/tabix/Makefile
- M /trunk/tabix/bgzip.c
- M /trunk/tabix/index.c
- M /trunk/tabix/main.c
-
- * tabix-0.0.0-1 (r500)
- * apparently working
-
-------------------------------------------------------------------------
-r499 | lh3lh3 | 2009-11-01 14:04:52 -0500 (Sun, 01 Nov 2009) | 2 lines
-Changed paths:
- D /trunk/tabix/parser.c
-
-obsolete file
-
-------------------------------------------------------------------------
-r498 | lh3lh3 | 2009-11-01 14:04:08 -0500 (Sun, 01 Nov 2009) | 2 lines
-Changed paths:
- M /trunk/tabix/bgzip.c
-
-bgzip is more like gzip in its command-line interface
-
-------------------------------------------------------------------------
-r497 | lh3lh3 | 2009-11-01 13:43:35 -0500 (Sun, 01 Nov 2009) | 2 lines
-Changed paths:
- A /trunk/tabix/Makefile
- A /trunk/tabix/bam_endian.h
- A /trunk/tabix/bgzf.c
- A /trunk/tabix/bgzf.h
- A /trunk/tabix/bgzip.c
- A /trunk/tabix/index.c
- A /trunk/tabix/khash.h
- A /trunk/tabix/knetfile.c
- A /trunk/tabix/knetfile.h
- A /trunk/tabix/ksort.h
- A /trunk/tabix/kstring.c
- A /trunk/tabix/kstring.h
- A /trunk/tabix/main.c
- A /trunk/tabix/parser.c
- A /trunk/tabix/tabix.h
-
-initial source code. It is BUGGY!
-
-------------------------------------------------------------------------
-r496 | lh3lh3 | 2009-11-01 13:42:39 -0500 (Sun, 01 Nov 2009) | 2 lines
-Changed paths:
- A /trunk/tabix
-
-A generic indexer for TAB-delimited genome position files
-
-------------------------------------------------------------------------
diff --git a/Makefile b/Makefile
index 0702ed3..7429ec0 100644
--- a/Makefile
+++ b/Makefile
@@ -1,71 +1,49 @@
-CC= gcc
-CPP= g++
-CFLAGS= -g -Wall -O2 -fPIC #-m64 #-arch ppc
+
+# Use ?= to allow override from the env or command-line.
+
+CC?= gcc
+CXX?= g++
+CXXFLAGS?= -g -Wall -O2 -fPIC #-m64 #-arch ppc
+INCLUDES?= -Ihtslib
+HTS_HEADERS?= htslib/htslib/bgzf.h htslib/htslib/tbx.h
+HTS_LIB?= htslib/libhts.a
+LIBPATH?= -L. -Lhtslib
+
DFLAGS= -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE
-LOBJS= bgzf.o kstring.o knetfile.o index.o bedidx.o
-AOBJS= main.o
-PROG= tabix bgzip tabix++
-INCLUDES=
-SUBDIRS= .
-LIBPATH=
-LIBCURSES=
+PROG= tabix++
+SUBDIRS=.
.SUFFIXES:.c .o
.c.o:
- $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $< -o $@
+ $(CC) -c $(CXXFLAGS) $(DFLAGS) $(INCLUDES) $< -o $@
all-recur lib-recur clean-recur cleanlocal-recur install-recur:
- @target=`echo $@ | sed s/-recur//`; \
- wdir=`pwd`; \
- list='$(SUBDIRS)'; for subdir in $$list; do \
- cd $$subdir; \
- $(MAKE) CC="$(CC)" DFLAGS="$(DFLAGS)" CFLAGS="$(CFLAGS)" \
- INCLUDES="$(INCLUDES)" LIBPATH="$(LIBPATH)" $$target || exit 1; \
- cd $$wdir; \
- done;
-
-all:$(PROG)
-
-lib:libtabix.a
-
-libtabix.so.1:$(LOBJS)
- $(CC) -shared -Wl,-soname,libtabix.so -o $@ $(LOBJS) -lc -lz
-
-libtabix.1.dylib:$(LOBJS)
- libtool -dynamic $(LOBJS) -o $@ -lc -lz
-
-libtabix.a:$(LOBJS)
- $(AR) -cru $@ $(LOBJS)
- ranlib $@
-
-tabix:lib $(AOBJS)
- $(CC) $(CFLAGS) -o $@ $(AOBJS) -lm $(LIBPATH) -L. -ltabix -lz
-
-tabix.o: tabix.hpp tabix.cpp
- $(CPP) $(CFLAGS) -c tabix.cpp
-
-tabix++:lib tabix.o main.cpp
- $(CPP) $(CFLAGS) -o $@ main.cpp tabix.o bgzf.o -lm $(LIBPATH) -L. -ltabix -lz
+ @target=`echo $@ | sed s/-recur//`; \
+ wdir=`pwd`; \
+ list='$(SUBDIRS)'; for subdir in $$list; do \
+ cd $$subdir; \
+ $(MAKE) CC="$(CC)" DFLAGS="$(DFLAGS)" CXXFLAGS="$(CXXFLAGS)" \
+ INCLUDES="$(INCLUDES)" LIBPATH="$(LIBPATH)" $$target \
+ || exit 1; \
+ cd $$wdir; \
+ done;
-bgzip:bgzip.o bgzf.o knetfile.o
- $(CC) $(CFLAGS) -o $@ bgzip.o bgzf.o knetfile.o -lz
+all: $(PROG)
-TabixReader.class:TabixReader.java
- javac -cp .:sam.jar TabixReader.java
+tabix.o: $(HTS_HEADERS) tabix.cpp tabix.hpp
+ $(CXX) $(CXXFLAGS) -c tabix.cpp $(INCLUDES)
-kstring.o:kstring.h
-knetfile.o:knetfile.h
-bgzf.o:bgzf.h knetfile.h
-index.o:bgzf.h tabix.h khash.h ksort.h kstring.h
-main.o:tabix.h kstring.h bgzf.h
-bgzip.o:bgzf.h
-bedidx.o:kseq.h khash.h
+htslib/libhts.a:
+ cd htslib && $(MAKE) lib-static
-tabix.pdf:tabix.tex
- pdflatex tabix.tex
+tabix++: tabix.o main.cpp $(HTS_LIB)
+ $(CXX) $(CXXFLAGS) -o $@ main.cpp tabix.o $(INCLUDES) $(LIBPATH) \
+ -lhts -lpthread -lm -lz
cleanlocal:
- rm -fr gmon.out *.o a.out *.dSYM $(PROG) *~ *.a tabix.aux tabix.log tabix.pdf *.class libtabix.*.dylib libtabix.so*
+ rm -fr gmon.out *.o a.out *.dSYM $(PROG) *~ *.a tabix.aux tabix.log \
+ tabix.pdf *.class libtabix.*.dylib libtabix.so*
+ cd htslib && $(MAKE) clean
clean:cleanlocal-recur
diff --git a/NEWS b/NEWS
deleted file mode 100644
index d230541..0000000
--- a/NEWS
+++ /dev/null
@@ -1,126 +0,0 @@
-Release 0.2.4 (10 April, 2011)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Notable changes:
-
- * Give an error if the index file is older than the data file.
-
- * Avoid a segfault given flawed input.
-
- * Added Python APIs contributed by Hyeshik Chang. The new APIs do not bind to
- the dynamic library and are reported to be faster. Pysam also comes with a
- tabix binding.
-
- * Added option "-r" for efficient header replacement.
-
- * Added BED support.
-
- * Synchronized the BGZF library between tabix and samtools.
-
-(0.2.4: 10 April 2011, r949)
-
-
-
-Beta Release 0.2.3 (8 December, 2010)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Notable changes:
-
- * Fixed a minor bug where the first record in a headerless file may be
- missed.
-
- * Added an option to print header lines.
-
- * Fixed a rare bug which may occasionally happen when retrieving data
- from a region without any records.
-
- * Enhanced error reporting.
-
- * Fixed a bug in bgzip which may delete the original file even if not
- intended.
-
-(0.2.3: 8 December 2010, r876)
-
-
-
-Beta Release 0.2.2 (28 June, 2010)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Notable changes:
-
- * Dropped the VCF3 support. Added VCF4 support.
-
- * Avoided the function name collision with samtools.
-
-(0.2.2: 28 June 2010, r603)
-
-
-
-Beta Release 0.2.1 (3 June, 2010)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Notable changes:
-
- * Allow shared library to be compiled. Added python binding to the
- shared library.
-
-(0.2.1: 3 June 2010, r582)
-
-
-
-Beta Release 0.2.0 (11 May, 2010)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Notable changes:
-
- * Fixed an issue for random access given an interval end larger than
- 2^29.
-
- * Updated the Java binding.
-
- * Added a Perl module using XS.
-
- * Improved the C APIs.
-
-(0.2.0: 11 May 2010, r574)
-
-
-
-Beta Release 0.1.6 (9 May, 2010)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Notable changes:
-
- * Improved backward compatibility. Release 0.1.5 does not work with the
- buggy index file generated by 0.1.2.
-
- * Fixed a bug in building linear index. The bug does not affect the
- results, only affects efficiency in rare cases.
-
- * Reduced the number of seek calls given an index generated by old
- version of tabix.
-
- * Added new APIs for retrieving data via an iterator. The old callback
- APIs are not changed, although internally it uses iterator to
- retrieve data.
-
-I am trying to freeze tabix. I just hope I am committing new bugs.
-
-(0.1.6: 9 May 2010, r563)
-
-
-
-Beta Release 0.1.5 (5 May, 2010)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Notable changes:
-
- * Clarified that tabix is released under MIT/X11.
-
- * Improved the robustness of indexing and retrieval.
-
- * Reduced the number of seek calls when the specified region starts
- from a 16kb block with no data. The index format is the same, but the
- content is changed a little.
-
-(0.1.5: 5 May 2010, r560)
diff --git a/README b/README
deleted file mode 100644
index 966660e..0000000
--- a/README
+++ /dev/null
@@ -1,6 +0,0 @@
-This is a fork of the tabix project [1] which includes a C++ class wrapper for
-reading tabix-indexed files.
-
-Author: Erik Garrison <erik.garrison at bc.edu>
-
-[1] http://samtools.sourceforge.net/tabix.shtml
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..aaa3f59
--- /dev/null
+++ b/README.md
@@ -0,0 +1,5 @@
+This is a C++ wrapper around [tabix
+project](http://samtools.sourceforge.net/tabix.shtml) which abstracts some
+of the details of opening and jumping in tabix-indexed files.
+
+Author: Erik Garrison <erik.garrison at gmail.com>
diff --git a/TabixReader.java b/TabixReader.java
deleted file mode 100644
index 5874202..0000000
--- a/TabixReader.java
+++ /dev/null
@@ -1,395 +0,0 @@
-/* The MIT License
-
- Copyright (c) 2010 Broad Institute.
-
- 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 <hengli at broadinstitute.org> */
-
-import net.sf.samtools.util.BlockCompressedInputStream;
-
-import java.io.*;
-import java.nio.*;
-import java.util.HashMap;
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.lang.StringBuffer;
-
-public class TabixReader
-{
- private String mFn;
- private BlockCompressedInputStream mFp;
-
- private int mPreset;
- private int mSc;
- private int mBc;
- private int mEc;
- private int mMeta;
- private int mSkip;
- private String[] mSeq;
-
- private HashMap<String, Integer> mChr2tid;
-
- private static int MAX_BIN = 37450;
- private static int TAD_MIN_CHUNK_GAP = 32768;
- private static int TAD_LIDX_SHIFT = 14;
-
- private class TPair64 implements Comparable<TPair64> {
- long u, v;
- public TPair64(final long _u, final long _v) {
- u = _u; v = _v;
- }
- public TPair64(final TPair64 p) {
- u = p.u; v = p.v;
- }
- public int compareTo(final TPair64 p) {
- return u == p.u? 0 : ((u < p.u) ^ (u < 0) ^ (p.u < 0))? -1 : 1; // unsigned 64-bit comparison
- }
- };
-
- private class TIndex {
- HashMap<Integer, TPair64[]> b; // binning index
- long[] l; // linear index
- };
- private TIndex[] mIndex;
-
- private class TIntv {
- int tid, beg, end;
- };
-
- private static boolean less64(final long u, final long v) { // unsigned 64-bit comparison
- return (u < v) ^ (u < 0) ^ (v < 0);
- }
-
- /**
- * The constructor
- *
- * @param fn File name of the data file
- */
- public TabixReader(final String fn) throws IOException {
- mFn = fn;
- mFp = new BlockCompressedInputStream(new File(fn));
- readIndex();
- }
-
- private static int reg2bins(final int beg, final int _end, final int[] list) {
- int i = 0, k, end = _end;
- if (beg >= end) return 0;
- if (end >= 1<<29) end = 1<<29;
- --end;
- list[i++] = 0;
- for (k = 1 + (beg>>26); k <= 1 + (end>>26); ++k) list[i++] = k;
- for (k = 9 + (beg>>23); k <= 9 + (end>>23); ++k) list[i++] = k;
- for (k = 73 + (beg>>20); k <= 73 + (end>>20); ++k) list[i++] = k;
- for (k = 585 + (beg>>17); k <= 585 + (end>>17); ++k) list[i++] = k;
- for (k = 4681 + (beg>>14); k <= 4681 + (end>>14); ++k) list[i++] = k;
- return i;
- }
-
- public static int readInt(final InputStream is) throws IOException {
- byte[] buf = new byte[4];
- is.read(buf);
- return ByteBuffer.wrap(buf).order(ByteOrder.LITTLE_ENDIAN).getInt();
- }
-
- public static long readLong(final InputStream is) throws IOException {
- byte[] buf = new byte[8];
- is.read(buf);
- return ByteBuffer.wrap(buf).order(ByteOrder.LITTLE_ENDIAN).getLong();
- }
-
- public static String readLine(final InputStream is) throws IOException {
- StringBuffer buf = new StringBuffer();
- int c;
- while ((c = is.read()) >= 0 && c != '\n')
- buf.append((char)c);
- if (c < 0) return null;
- return buf.toString();
- }
-
- /**
- * Read the Tabix index from a file
- *
- * @param fp File pointer
- */
- public void readIndex(final File fp) throws IOException {
- if (fp == null) return;
- BlockCompressedInputStream is = new BlockCompressedInputStream(fp);
- byte[] buf = new byte[4];
-
- is.read(buf, 0, 4); // read "TBI\1"
- mSeq = new String[readInt(is)]; // # sequences
- mChr2tid = new HashMap<String, Integer>();
- mPreset = readInt(is);
- mSc = readInt(is);
- mBc = readInt(is);
- mEc = readInt(is);
- mMeta = readInt(is);
- mSkip = readInt(is);
- // read sequence dictionary
- int i, j, k, l = readInt(is);
- buf = new byte[l];
- is.read(buf);
- for (i = j = k = 0; i < buf.length; ++i) {
- if (buf[i] == 0) {
- byte[] b = new byte[i - j];
- System.arraycopy(buf, j, b, 0, b.length);
- String s = new String(b);
- mChr2tid.put(s, k);
- mSeq[k++] = s;
- j = i + 1;
- }
- }
- // read the index
- mIndex = new TIndex[mSeq.length];
- for (i = 0; i < mSeq.length; ++i) {
- // the binning index
- int n_bin = readInt(is);
- mIndex[i] = new TIndex();
- mIndex[i].b = new HashMap<Integer, TPair64[]>();
- for (j = 0; j < n_bin; ++j) {
- int bin = readInt(is);
- TPair64[] chunks = new TPair64[readInt(is)];
- for (k = 0; k < chunks.length; ++k) {
- long u = readLong(is);
- long v = readLong(is);
- chunks[k] = new TPair64(u, v); // in C, this is inefficient
- }
- mIndex[i].b.put(bin, chunks);
- }
- // the linear index
- mIndex[i].l = new long[readInt(is)];
- for (k = 0; k < mIndex[i].l.length; ++k)
- mIndex[i].l[k] = readLong(is);
- }
- // close
- is.close();
- }
-
- /**
- * Read the Tabix index from the default file.
- */
- public void readIndex() throws IOException {
- readIndex(new File(mFn + ".tbi"));
- }
-
- /**
- * Read one line from the data file.
- */
- public String readLine() throws IOException {
- return readLine(mFp);
- }
-
- private int chr2tid(final String chr) {
- if (mChr2tid.containsKey(chr)) return mChr2tid.get(chr);
- else return -1;
- }
-
- /**
- * Parse a region in the format of "chr1", "chr1:100" or "chr1:100-1000"
- *
- * @param reg Region string
- * @return An array where the three elements are sequence_id,
- * region_begin and region_end. On failure, sequence_id==-1.
- */
- public int[] parseReg(final String reg) { // FIXME: NOT working when the sequence name contains : or -.
- String chr;
- int colon, hyphen;
- int[] ret = new int[3];
- colon = reg.indexOf(':'); hyphen = reg.indexOf('-');
- chr = colon >= 0? reg.substring(0, colon) : reg;
- ret[1] = colon >= 0? Integer.parseInt(reg.substring(colon+1, hyphen >= 0? hyphen : reg.length())) - 1 : 0;
- ret[2] = hyphen >= 0? Integer.parseInt(reg.substring(hyphen+1)) : 0x7fffffff;
- ret[0] = chr2tid(chr);
- return ret;
- }
-
- private TIntv getIntv(final String s) {
- TIntv intv = new TIntv();
- int col = 0, end = 0, beg = 0;
- while ((end = s.indexOf('\t', beg)) >= 0 || end == -1) {
- ++col;
- if (col == mSc) {
- intv.tid = chr2tid(s.substring(beg, end));
- } else if (col == mBc) {
- intv.beg = intv.end = Integer.parseInt(s.substring(beg, end));
- if ((mPreset&0x10000) != 0) ++intv.end;
- else --intv.beg;
- if (intv.beg < 0) intv.beg = 0;
- if (intv.end < 1) intv.end = 1;
- } else { // FIXME: SAM supports are not tested yet
- if ((mPreset&0xffff) == 0) { // generic
- if (col == mEc)
- intv.end = Integer.parseInt(s.substring(beg, end));
- } else if ((mPreset&0xffff) == 1) { // SAM
- if (col == 6) { // CIGAR
- int l = 0, i, j;
- String cigar = s.substring(beg, end);
- for (i = j = 0; i < cigar.length(); ++i) {
- if (cigar.charAt(i) > '9') {
- int op = cigar.charAt(i);
- if (op == 'M' || op == 'D' || op == 'N')
- l += Integer.parseInt(cigar.substring(j, i));
- }
- }
- intv.end = intv.beg + l;
- }
- } else if ((mPreset&0xffff) == 2) { // VCF
- String alt;
- alt = end >= 0? s.substring(beg, end) : s.substring(beg);
- if (col == 4) { // REF
- if (alt.length() > 0) intv.end = intv.beg + alt.length();
- } else if (col == 8) { // INFO
- int e_off = -1, i = alt.indexOf("END=");
- if (i == 0) e_off = 4;
- else if (i > 0) {
- i = alt.indexOf(";END=");
- if (i >= 0) e_off = i + 5;
- }
- if (e_off > 0) {
- i = alt.indexOf(";", e_off);
- intv.end = Integer.parseInt(i > e_off? alt.substring(e_off, i) : alt.substring(e_off));
- }
- }
- }
- }
- if (end == -1) break;
- beg = end + 1;
- }
- return intv;
- }
-
- public class Iterator {
- private int i, n_seeks;
- private int tid, beg, end;
- private TPair64[] off;
- private long curr_off;
- private boolean iseof;
-
- public Iterator(final int _tid, final int _beg, final int _end, final TPair64[] _off) {
- i = -1; n_seeks = 0; curr_off = 0; iseof = false;
- off = _off; tid = _tid; beg = _beg; end = _end;
- }
-
- public String next() throws IOException {
- if (iseof) return null;
- for (;;) {
- if (curr_off == 0 || !less64(curr_off, off[i].v)) { // then jump to the next chunk
- if (i == off.length - 1) break; // no more chunks
- if (i >= 0) assert(curr_off == off[i].v); // otherwise bug
- if (i < 0 || off[i].v != off[i+1].u) { // not adjacent chunks; then seek
- mFp.seek(off[i+1].u);
- curr_off = mFp.getFilePointer();
- ++n_seeks;
- }
- ++i;
- }
- String s;
- if ((s = readLine(mFp)) != null) {
- TIntv intv;
- char[] str = s.toCharArray();
- curr_off = mFp.getFilePointer();
- if (str.length == 0 || str[0] == mMeta) continue;
- intv = getIntv(s);
- if (intv.tid != tid || intv.beg >= end) break; // no need to proceed
- else if (intv.end > beg && intv.beg < end) return s; // overlap; return
- } else break; // end of file
- }
- iseof = true;
- return null;
- }
- };
-
- public Iterator query(final int tid, final int beg, final int end) {
- TPair64[] off, chunks;
- long min_off;
- TIndex idx = mIndex[tid];
- int[] bins = new int[MAX_BIN];
- int i, l, n_off, n_bins = reg2bins(beg, end, bins);
- if (idx.l.length > 0)
- min_off = (beg>>TAD_LIDX_SHIFT >= idx.l.length)? idx.l[idx.l.length-1] : idx.l[beg>>TAD_LIDX_SHIFT];
- else min_off = 0;
- for (i = n_off = 0; i < n_bins; ++i) {
- if ((chunks = idx.b.get(bins[i])) != null)
- n_off += chunks.length;
- }
- if (n_off == 0) return null;
- off = new TPair64[n_off];
- for (i = n_off = 0; i < n_bins; ++i)
- if ((chunks = idx.b.get(bins[i])) != null)
- for (int j = 0; j < chunks.length; ++j)
- if (less64(min_off, chunks[j].v))
- off[n_off++] = new TPair64(chunks[j]);
- if (n_off == 0) return null;
- Arrays.sort(off, 0, n_off);
- // resolve completely contained adjacent blocks
- for (i = 1, l = 0; i < n_off; ++i) {
- if (less64(off[l].v, off[i].v)) {
- ++l;
- off[l].u = off[i].u; off[l].v = off[i].v;
- }
- }
- n_off = l + 1;
- // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
- for (i = 1; i < n_off; ++i)
- if (!less64(off[i-1].v, off[i].u)) off[i-1].v = off[i].u;
- // merge adjacent blocks
- for (i = 1, l = 0; i < n_off; ++i) {
- if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v;
- else {
- ++l;
- off[l].u = off[i].u;
- off[l].v = off[i].v;
- }
- }
- n_off = l + 1;
- // return
- TPair64[] ret = new TPair64[n_off];
- for (i = 0; i < n_off; ++i) ret[i] = new TPair64(off[i].u, off[i].v); // in C, this is inefficient
- return new TabixReader.Iterator(tid, beg, end, ret);
- }
-
- public Iterator query(final String reg) {
- int[] x = parseReg(reg);
- return query(x[0], x[1], x[2]);
- }
-
- public static void main(String[] args) {
- if (args.length < 1) {
- System.out.println("Usage: java -cp .:sam.jar TabixReader <in.gz> [region]");
- System.exit(1);
- }
- try {
- TabixReader tr = new TabixReader(args[0]);
- String s;
- if (args.length == 1) { // no region is specified; print the whole file
- while ((s = tr.readLine()) != null)
- System.out.println(s);
- } else { // a region is specified; random access
- TabixReader.Iterator iter = tr.query(args[1]); // get the iterator
- while (iter != null && (s = iter.next()) != null)
- System.out.println(s);
- }
- } catch (IOException e) {
- }
- }
-}
diff --git a/bam_endian.h b/bam_endian.h
deleted file mode 100644
index 0fc74a8..0000000
--- a/bam_endian.h
+++ /dev/null
@@ -1,42 +0,0 @@
-#ifndef BAM_ENDIAN_H
-#define BAM_ENDIAN_H
-
-#include <stdint.h>
-
-static inline int bam_is_big_endian()
-{
- long one= 1;
- return !(*((char *)(&one)));
-}
-static inline uint16_t bam_swap_endian_2(uint16_t v)
-{
- return (uint16_t)(((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8));
-}
-static inline void *bam_swap_endian_2p(void *x)
-{
- *(uint16_t*)x = bam_swap_endian_2(*(uint16_t*)x);
- return x;
-}
-static inline uint32_t bam_swap_endian_4(uint32_t v)
-{
- v = ((v & 0x0000FFFFU) << 16) | (v >> 16);
- return ((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8);
-}
-static inline void *bam_swap_endian_4p(void *x)
-{
- *(uint32_t*)x = bam_swap_endian_4(*(uint32_t*)x);
- return x;
-}
-static inline uint64_t bam_swap_endian_8(uint64_t v)
-{
- v = ((v & 0x00000000FFFFFFFFLLU) << 32) | (v >> 32);
- v = ((v & 0x0000FFFF0000FFFFLLU) << 16) | ((v & 0xFFFF0000FFFF0000LLU) >> 16);
- return ((v & 0x00FF00FF00FF00FFLLU) << 8) | ((v & 0xFF00FF00FF00FF00LLU) >> 8);
-}
-static inline void *bam_swap_endian_8p(void *x)
-{
- *(uint64_t*)x = bam_swap_endian_8(*(uint64_t*)x);
- return x;
-}
-
-#endif
diff --git a/bedidx.c b/bedidx.c
deleted file mode 100644
index 722877d..0000000
--- a/bedidx.c
+++ /dev/null
@@ -1,156 +0,0 @@
-#include <stdlib.h>
-#include <stdint.h>
-#include <string.h>
-#include <stdio.h>
-#include <zlib.h>
-
-#include "ksort.h"
-KSORT_INIT_GENERIC(uint64_t)
-
-#include "kseq.h"
-KSTREAM_INIT(gzFile, gzread, 8192)
-
-typedef struct {
- int n, m;
- uint64_t *a;
- int *idx;
-} bed_reglist_t;
-
-#include "khash.h"
-KHASH_MAP_INIT_STR(reg, bed_reglist_t)
-
-#define LIDX_SHIFT 13
-
-typedef kh_reg_t reghash_t;
-
-int *bed_index_core(int n, uint64_t *a, int *n_idx)
-{
- int i, j, m, *idx;
- m = *n_idx = 0; idx = 0;
- for (i = 0; i < n; ++i) {
- int beg, end;
- beg = a[i]>>32 >> LIDX_SHIFT; end = ((uint32_t)a[i]) >> LIDX_SHIFT;
- if (m < end + 1) {
- int oldm = m;
- m = end + 1;
- kroundup32(m);
- idx = realloc(idx, m * sizeof(int));
- for (j = oldm; j < m; ++j) idx[j] = -1;
- }
- if (beg == end) {
- if (idx[beg] < 0) idx[beg] = i;
- } else {
- for (j = beg; j <= end; ++j)
- if (idx[j] < 0) idx[j] = i;
- }
- *n_idx = end + 1;
- }
- return idx;
-}
-
-void bed_index(void *_h)
-{
- reghash_t *h = (reghash_t*)_h;
- khint_t k;
- for (k = 0; k < kh_end(h); ++k) {
- if (kh_exist(h, k)) {
- bed_reglist_t *p = &kh_val(h, k);
- if (p->idx) free(p->idx);
- ks_introsort(uint64_t, p->n, p->a);
- p->idx = bed_index_core(p->n, p->a, &p->m);
- }
- }
-}
-
-int bed_overlap_core(const bed_reglist_t *p, int beg, int end)
-{
- int i, min_off;
- if (p->n == 0) return 0;
- min_off = (beg>>LIDX_SHIFT >= p->n)? p->idx[p->n-1] : p->idx[beg>>LIDX_SHIFT];
- if (min_off < 0) { // TODO: this block can be improved, but speed should not matter too much here
- int n = beg>>LIDX_SHIFT;
- if (n > p->n) n = p->n;
- for (i = n - 1; i >= 0; --i)
- if (p->idx[i] >= 0) break;
- min_off = i >= 0? p->idx[i] : 0;
- }
- for (i = min_off; i < p->n; ++i) {
- if ((int)(p->a[i]>>32) >= end) break; // out of range; no need to proceed
- if ((int32_t)p->a[i] > beg && (int32_t)(p->a[i]>>32) < end)
- return 1; // find the overlap; return
- }
- return 0;
-}
-
-int bed_overlap(const void *_h, const char *chr, int beg, int end)
-{
- const reghash_t *h = (const reghash_t*)_h;
- khint_t k;
- if (!h) return 0;
- k = kh_get(reg, h, chr);
- if (k == kh_end(h)) return 0;
- return bed_overlap_core(&kh_val(h, k), beg, end);
-}
-
-void *bed_read(const char *fn)
-{
- reghash_t *h = kh_init(reg);
- gzFile fp;
- kstream_t *ks;
- int dret;
- kstring_t *str;
- // read the list
- fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
- if (fp == 0) return 0;
- str = calloc(1, sizeof(kstring_t));
- ks = ks_init(fp);
- while (ks_getuntil(ks, 0, str, &dret) >= 0) { // read the chr name
- int beg = -1, end = -1;
- bed_reglist_t *p;
- khint_t k = kh_get(reg, h, str->s);
- if (k == kh_end(h)) { // absent from the hash table
- int ret;
- char *s = strdup(str->s);
- k = kh_put(reg, h, s, &ret);
- memset(&kh_val(h, k), 0, sizeof(bed_reglist_t));
- }
- p = &kh_val(h, k);
- if (dret != '\n') { // if the lines has other characters
- if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
- beg = atoi(str->s); // begin
- if (dret != '\n') {
- if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0]))
- end = atoi(str->s); // end
- }
- }
- }
- if (dret != '\n') while ((dret = ks_getc(ks)) > 0 && dret != '\n'); // skip the rest of the line
- if (end < 0 && beg > 0) end = beg, beg = beg - 1; // if there is only one column
- if (beg >= 0 && end > beg) {
- if (p->n == p->m) {
- p->m = p->m? p->m<<1 : 4;
- p->a = realloc(p->a, p->m * 8);
- }
- p->a[p->n++] = (uint64_t)beg<<32 | end;
- }
- }
- ks_destroy(ks);
- gzclose(fp);
- free(str->s); free(str);
- bed_index(h);
- return h;
-}
-
-void bed_destroy(void *_h)
-{
- reghash_t *h = (reghash_t*)_h;
- khint_t k;
- for (k = 0; k < kh_end(h); ++k) {
- if (kh_exist(h, k)) {
- free(kh_val(h, k).a);
- free(kh_val(h, k).idx);
- free((char*)kh_key(h, k));
- }
- }
- kh_destroy(reg, h);
-}
diff --git a/bgzf.c b/bgzf.c
deleted file mode 100644
index 216cd04..0000000
--- a/bgzf.c
+++ /dev/null
@@ -1,714 +0,0 @@
-/* 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-25 by lh3: optionally use my knetfile library to access file on a FTP.
- 2009-06-12 by lh3: support a mode string like "wu" where 'u' for uncompressed output */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <unistd.h>
-#include <fcntl.h>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include "bgzf.h"
-
-#include "khash.h"
-typedef struct {
- int size;
- uint8_t *block;
- int64_t end_offset;
-} 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);
-#endif
-
-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 BLOCK_HEADER_LENGTH = 18;
-static const int BLOCK_FOOTER_LENGTH = 8;
-
-static const int GZIP_ID1 = 31;
-static const int GZIP_ID2 = 139;
-static const int CM_DEFLATE = 8;
-static const int FLG_FEXTRA = 4;
-static const int OS_UNKNOWN = 255;
-static const int BGZF_ID1 = 66; // 'B'
-static const int BGZF_ID2 = 67; // 'C'
-static const int BGZF_LEN = 2;
-static const int BGZF_XLEN = 6; // BGZF_LEN+4
-
-static const int GZIP_WINDOW_BITS = -15; // no zlib header
-static const int Z_DEFAULT_MEM_LEVEL = 8;
-
-
-inline
-void
-packInt16(uint8_t* buffer, uint16_t value)
-{
- buffer[0] = value;
- buffer[1] = value >> 8;
-}
-
-inline
-int
-unpackInt16(const uint8_t* buffer)
-{
- return (buffer[0] | (buffer[1] << 8));
-}
-
-inline
-void
-packInt32(uint8_t* buffer, uint32_t value)
-{
- buffer[0] = value;
- buffer[1] = value >> 8;
- buffer[2] = value >> 16;
- buffer[3] = value >> 24;
-}
-
-static inline
-int
-bgzf_min(int x, int y)
-{
- return (x < y) ? x : y;
-}
-
-static
-void
-report_error(BGZF* fp, const char* message) {
- fp->error = message;
-}
-
-int bgzf_check_bgzf(const char *fn)
-{
- BGZF *fp;
- uint8_t buf[10],magic[10]="\037\213\010\4\0\0\0\0\0\377";
- int n;
-
- if ((fp = bgzf_open(fn, "r")) == 0)
- {
- fprintf(stderr, "[bgzf_check_bgzf] failed to open the file: %s\n",fn);
- return -1;
- }
-
-#ifdef _USE_KNETFILE
- n = knet_read(fp->x.fpr, buf, 10);
-#else
- n = fread(buf, 1, 10, fp->file);
-#endif
- bgzf_close(fp);
-
- if ( n!=10 )
- return -1;
-
- if ( !memcmp(magic, buf, 10) ) return 1;
- return 0;
-}
-
-static BGZF *bgzf_read_init()
-{
- BGZF *fp;
- fp = calloc(1, sizeof(BGZF));
- fp->uncompressed_block_size = MAX_BLOCK_SIZE;
- fp->uncompressed_block = malloc(MAX_BLOCK_SIZE);
- fp->compressed_block_size = MAX_BLOCK_SIZE;
- fp->compressed_block = malloc(MAX_BLOCK_SIZE);
- fp->cache_size = 0;
- fp->cache = kh_init(cache);
- return fp;
-}
-
-static
-BGZF*
-open_read(int fd)
-{
-#ifdef _USE_KNETFILE
- knetFile *file = knet_dopen(fd, "r");
-#else
- FILE* file = fdopen(fd, "r");
-#endif
- BGZF* fp;
- if (file == 0) return 0;
- fp = bgzf_read_init();
- fp->file_descriptor = fd;
- fp->open_mode = 'r';
-#ifdef _USE_KNETFILE
- fp->x.fpr = file;
-#else
- fp->file = file;
-#endif
- return fp;
-}
-
-static
-BGZF*
-open_write(int fd, int compress_level) // compress_level==-1 for the default level
-{
- FILE* file = fdopen(fd, "w");
- BGZF* fp;
- if (file == 0) return 0;
- fp = malloc(sizeof(BGZF));
- fp->file_descriptor = fd;
- fp->open_mode = 'w';
- fp->owned_file = 0;
- fp->compress_level = compress_level < 0? Z_DEFAULT_COMPRESSION : compress_level; // Z_DEFAULT_COMPRESSION==-1
- if (fp->compress_level > 9) fp->compress_level = Z_DEFAULT_COMPRESSION;
-#ifdef _USE_KNETFILE
- fp->x.fpw = file;
-#else
- fp->file = file;
-#endif
- fp->uncompressed_block_size = DEFAULT_BLOCK_SIZE;
- fp->uncompressed_block = NULL;
- fp->compressed_block_size = MAX_BLOCK_SIZE;
- fp->compressed_block = malloc(MAX_BLOCK_SIZE);
- fp->block_address = 0;
- fp->block_offset = 0;
- fp->block_length = 0;
- fp->error = NULL;
- return fp;
-}
-
-BGZF*
-bgzf_open(const char* __restrict path, const char* __restrict mode)
-{
- BGZF* fp = NULL;
- if (strchr(mode, 'r') || strchr(mode, 'R')) { /* The reading mode is preferred. */
-#ifdef _USE_KNETFILE
- knetFile *file = knet_open(path, mode);
- if (file == 0) return 0;
- fp = bgzf_read_init();
- fp->file_descriptor = -1;
- fp->open_mode = 'r';
- fp->x.fpr = file;
-#else
- 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 (strchr(mode, 'w') || strchr(mode, 'W')) {
- int fd, compress_level = -1, oflag = O_WRONLY | O_CREAT | O_TRUNC;
-#ifdef _WIN32
- oflag |= O_BINARY;
-#endif
- fd = open(path, oflag, 0666);
- if (fd == -1) return 0;
- { // set compress_level
- int i;
- for (i = 0; mode[i]; ++i)
- if (mode[i] >= '0' && mode[i] <= '9') break;
- if (mode[i]) compress_level = (int)mode[i] - '0';
- if (strchr(mode, 'u')) compress_level = 0;
- }
- fp = open_write(fd, compress_level);
- }
- if (fp != NULL) fp->owned_file = 1;
- return fp;
-}
-
-BGZF*
-bgzf_fdopen(int fd, const char * __restrict mode)
-{
- if (fd == -1) return 0;
- if (mode[0] == 'r' || mode[0] == 'R') {
- return open_read(fd);
- } else if (mode[0] == 'w' || mode[0] == 'W') {
- int i, compress_level = -1;
- for (i = 0; mode[i]; ++i)
- if (mode[i] >= '0' && mode[i] <= '9') break;
- if (mode[i]) compress_level = (int)mode[i] - '0';
- if (strchr(mode, 'u')) compress_level = 0;
- return open_write(fd, compress_level);
- } else {
- return NULL;
- }
-}
-
-static
-int
-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.
-
- bgzf_byte_t* buffer = fp->compressed_block;
- int buffer_size = fp->compressed_block_size;
-
- // Init gzip header
- buffer[0] = GZIP_ID1;
- buffer[1] = GZIP_ID2;
- buffer[2] = CM_DEFLATE;
- buffer[3] = FLG_FEXTRA;
- buffer[4] = 0; // mtime
- buffer[5] = 0;
- buffer[6] = 0;
- buffer[7] = 0;
- buffer[8] = 0;
- buffer[9] = OS_UNKNOWN;
- buffer[10] = BGZF_XLEN;
- buffer[11] = 0;
- buffer[12] = BGZF_ID1;
- buffer[13] = BGZF_ID2;
- buffer[14] = BGZF_LEN;
- buffer[15] = 0;
- buffer[16] = 0; // placeholder for block length
- buffer[17] = 0;
-
- // loop to retry for blocks that do not compress enough
- int input_length = block_length;
- int compressed_length = 0;
- while (1) {
- z_stream zs;
- zs.zalloc = NULL;
- zs.zfree = NULL;
- zs.next_in = fp->uncompressed_block;
- zs.avail_in = input_length;
- zs.next_out = (void*)&buffer[BLOCK_HEADER_LENGTH];
- zs.avail_out = buffer_size - BLOCK_HEADER_LENGTH - BLOCK_FOOTER_LENGTH;
-
- int status = deflateInit2(&zs, fp->compress_level, Z_DEFLATED,
- GZIP_WINDOW_BITS, Z_DEFAULT_MEM_LEVEL, Z_DEFAULT_STRATEGY);
- if (status != Z_OK) {
- report_error(fp, "deflate init failed");
- return -1;
- }
- status = deflate(&zs, Z_FINISH);
- if (status != Z_STREAM_END) {
- deflateEnd(&zs);
- if (status == Z_OK) {
- // Not enough space in buffer.
- // Can happen in the rare case the input doesn't compress enough.
- // Reduce the amount of input until it fits.
- input_length -= 1024;
- if (input_length <= 0) {
- // should never happen
- report_error(fp, "input reduction failed");
- return -1;
- }
- continue;
- }
- report_error(fp, "deflate failed");
- return -1;
- }
- status = deflateEnd(&zs);
- if (status != Z_OK) {
- report_error(fp, "deflate end failed");
- return -1;
- }
- compressed_length = zs.total_out;
- compressed_length += BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;
- if (compressed_length > MAX_BLOCK_SIZE) {
- // should never happen
- report_error(fp, "deflate overflow");
- return -1;
- }
- break;
- }
-
- packInt16((uint8_t*)&buffer[16], compressed_length-1);
- uint32_t crc = crc32(0L, NULL, 0L);
- crc = crc32(crc, fp->uncompressed_block, input_length);
- packInt32((uint8_t*)&buffer[compressed_length-8], crc);
- packInt32((uint8_t*)&buffer[compressed_length-4], input_length);
-
- int remaining = block_length - input_length;
- if (remaining > 0) {
- if (remaining > input_length) {
- // should never happen (check so we can use memcpy)
- report_error(fp, "remainder too large");
- return -1;
- }
- memcpy(fp->uncompressed_block,
- fp->uncompressed_block + input_length,
- remaining);
- }
- fp->block_offset = remaining;
- return compressed_length;
-}
-
-static
-int
-inflate_block(BGZF* fp, int block_length)
-{
- // Inflate the block in fp->compressed_block into fp->uncompressed_block
-
- z_stream zs;
- int status;
- zs.zalloc = NULL;
- zs.zfree = NULL;
- zs.next_in = fp->compressed_block + 18;
- zs.avail_in = block_length - 16;
- zs.next_out = fp->uncompressed_block;
- zs.avail_out = fp->uncompressed_block_size;
-
- status = inflateInit2(&zs, GZIP_WINDOW_BITS);
- if (status != Z_OK) {
- report_error(fp, "inflate init failed");
- return -1;
- }
- status = inflate(&zs, Z_FINISH);
- if (status != Z_STREAM_END) {
- inflateEnd(&zs);
- report_error(fp, "inflate failed");
- return -1;
- }
- status = inflateEnd(&zs);
- if (status != Z_OK) {
- report_error(fp, "inflate failed");
- return -1;
- }
- return zs.total_out;
-}
-
-static
-int
-check_header(const bgzf_byte_t* header)
-{
- return (header[0] == GZIP_ID1 &&
- header[1] == (bgzf_byte_t) GZIP_ID2 &&
- header[2] == Z_DEFLATED &&
- (header[3] & FLG_FEXTRA) != 0 &&
- unpackInt16((uint8_t*)&header[10]) == BGZF_XLEN &&
- header[12] == BGZF_ID1 &&
- header[13] == BGZF_ID2 &&
- unpackInt16((uint8_t*)&header[14]) == BGZF_LEN);
-}
-
-static void free_cache(BGZF *fp)
-{
- khint_t k;
- khash_t(cache) *h = (khash_t(cache)*)fp->cache;
- if (fp->open_mode != 'r') return;
- for (k = kh_begin(h); k < kh_end(h); ++k)
- if (kh_exist(h, k)) free(kh_val(h, k).block);
- kh_destroy(cache, h);
-}
-
-static int load_block_from_cache(BGZF *fp, int64_t block_address)
-{
- khint_t k;
- cache_t *p;
- khash_t(cache) *h = (khash_t(cache)*)fp->cache;
- k = kh_get(cache, h, block_address);
- if (k == kh_end(h)) return 0;
- p = &kh_val(h, k);
- if (fp->block_length != 0) fp->block_offset = 0;
- fp->block_address = block_address;
- fp->block_length = p->size;
- memcpy(fp->uncompressed_block, p->block, MAX_BLOCK_SIZE);
-#ifdef _USE_KNETFILE
- knet_seek(fp->x.fpr, p->end_offset, SEEK_SET);
-#else
- fseeko(fp->file, p->end_offset, SEEK_SET);
-#endif
- return p->size;
-}
-
-static void cache_block(BGZF *fp, int size)
-{
- int ret;
- khint_t k;
- cache_t *p;
- khash_t(cache) *h = (khash_t(cache)*)fp->cache;
- if (MAX_BLOCK_SIZE >= fp->cache_size) return;
- if ((kh_size(h) + 1) * MAX_BLOCK_SIZE > fp->cache_size) {
- /* A better way would be to remove the oldest block in the
- * cache, but here we remove a random one for simplicity. This
- * should not have a big impact on performance. */
- for (k = kh_begin(h); k < kh_end(h); ++k)
- if (kh_exist(h, k)) break;
- if (k < kh_end(h)) {
- free(kh_val(h, k).block);
- kh_del(cache, h, k);
- }
- }
- k = kh_put(cache, h, fp->block_address, &ret);
- if (ret == 0) return; // if this happens, a bug!
- p = &kh_val(h, k);
- p->size = fp->block_length;
- p->end_offset = fp->block_address + size;
- p->block = malloc(MAX_BLOCK_SIZE);
- memcpy(kh_val(h, k).block, fp->uncompressed_block, MAX_BLOCK_SIZE);
-}
-
-int
-bgzf_read_block(BGZF* fp)
-{
- bgzf_byte_t header[BLOCK_HEADER_LENGTH];
- int count, size = 0, block_length, remaining;
-#ifdef _USE_KNETFILE
- int64_t block_address = knet_tell(fp->x.fpr);
- if (load_block_from_cache(fp, block_address)) return 0;
- count = knet_read(fp->x.fpr, header, sizeof(header));
-#else
- int64_t block_address = ftello(fp->file);
- if (load_block_from_cache(fp, block_address)) return 0;
- count = fread(header, 1, sizeof(header), fp->file);
-#endif
- if (count == 0) {
- fp->block_length = 0;
- return 0;
- }
- size = count;
- if (count != sizeof(header)) {
- report_error(fp, "read failed");
- return -1;
- }
- if (!check_header(header)) {
- report_error(fp, "invalid block header");
- return -1;
- }
- block_length = unpackInt16((uint8_t*)&header[16]) + 1;
- bgzf_byte_t* compressed_block = (bgzf_byte_t*) fp->compressed_block;
- memcpy(compressed_block, header, BLOCK_HEADER_LENGTH);
- remaining = block_length - BLOCK_HEADER_LENGTH;
-#ifdef _USE_KNETFILE
- count = knet_read(fp->x.fpr, &compressed_block[BLOCK_HEADER_LENGTH], remaining);
-#else
- count = fread(&compressed_block[BLOCK_HEADER_LENGTH], 1, remaining, fp->file);
-#endif
- if (count != remaining) {
- report_error(fp, "read failed");
- return -1;
- }
- size += count;
- count = inflate_block(fp, block_length);
- if (count < 0) return -1;
- if (fp->block_length != 0) {
- // Do not reset offset if this read follows a seek.
- fp->block_offset = 0;
- }
- fp->block_address = block_address;
- fp->block_length = count;
- cache_block(fp, size);
- return 0;
-}
-
-int
-bgzf_read(BGZF* fp, void* data, int length)
-{
- if (length <= 0) {
- return 0;
- }
- if (fp->open_mode != 'r') {
- report_error(fp, "file not open for reading");
- return -1;
- }
-
- int bytes_read = 0;
- bgzf_byte_t* output = data;
- while (bytes_read < length) {
- int copy_length, available = fp->block_length - fp->block_offset;
- bgzf_byte_t *buffer;
- if (available <= 0) {
- if (bgzf_read_block(fp) != 0) {
- return -1;
- }
- available = fp->block_length - fp->block_offset;
- if (available <= 0) {
- break;
- }
- }
- copy_length = bgzf_min(length-bytes_read, available);
- buffer = fp->uncompressed_block;
- memcpy(output, buffer + fp->block_offset, copy_length);
- fp->block_offset += copy_length;
- output += copy_length;
- bytes_read += copy_length;
- }
- if (fp->block_offset == fp->block_length) {
-#ifdef _USE_KNETFILE
- fp->block_address = knet_tell(fp->x.fpr);
-#else
- fp->block_address = ftello(fp->file);
-#endif
- fp->block_offset = 0;
- fp->block_length = 0;
- }
- return bytes_read;
-}
-
-int bgzf_flush(BGZF* fp)
-{
- while (fp->block_offset > 0) {
- int count, block_length;
- block_length = deflate_block(fp, fp->block_offset);
- if (block_length < 0) return -1;
-#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
- if (count != block_length) {
- report_error(fp, "write failed");
- return -1;
- }
- fp->block_address += block_length;
- }
- return 0;
-}
-
-int bgzf_flush_try(BGZF *fp, int size)
-{
- if (fp->block_offset + size > fp->uncompressed_block_size)
- return bgzf_flush(fp);
- return -1;
-}
-
-int bgzf_write(BGZF* fp, const void* data, int length)
-{
- const bgzf_byte_t *input = data;
- int block_length, bytes_written;
- if (fp->open_mode != 'w') {
- report_error(fp, "file not open for writing");
- return -1;
- }
-
- if (fp->uncompressed_block == NULL)
- fp->uncompressed_block = malloc(fp->uncompressed_block_size);
-
- input = data;
- block_length = fp->uncompressed_block_size;
- bytes_written = 0;
- while (bytes_written < length) {
- 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;
- bytes_written += copy_length;
- if (fp->block_offset == block_length) {
- if (bgzf_flush(fp) != 0) {
- break;
- }
- }
- }
- return bytes_written;
-}
-
-int bgzf_close(BGZF* fp)
-{
- if (fp->open_mode == 'w') {
- if (bgzf_flush(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
- if (fflush(fp->file) != 0) {
-#endif
- report_error(fp, "flush failed");
- return -1;
- }
- }
- if (fp->owned_file) {
-#ifdef _USE_KNETFILE
- int ret;
- if (fp->open_mode == 'w') ret = fclose(fp->x.fpw);
- else ret = knet_close(fp->x.fpr);
- if (ret != 0) return -1;
-#else
- if (fclose(fp->file) != 0) return -1;
-#endif
- }
- free(fp->uncompressed_block);
- free(fp->compressed_block);
- free_cache(fp);
- free(fp);
- return 0;
-}
-
-void bgzf_set_cache_size(BGZF *fp, int 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)
-{
- int block_offset;
- int64_t block_address;
-
- if (fp->open_mode != 'r') {
- report_error(fp, "file not open for read");
- return -1;
- }
- if (where != SEEK_SET) {
- report_error(fp, "unimplemented seek option");
- return -1;
- }
- block_offset = pos & 0xFFFF;
- block_address = (pos >> 16) & 0xFFFFFFFFFFFFLL;
-#ifdef _USE_KNETFILE
- if (knet_seek(fp->x.fpr, block_address, SEEK_SET) != 0) {
-#else
- if (fseeko(fp->file, block_address, SEEK_SET) != 0) {
-#endif
- report_error(fp, "seek failed");
- return -1;
- }
- fp->block_length = 0; // indicates current block is not loaded
- fp->block_address = block_address;
- fp->block_offset = block_offset;
- return 0;
-}
diff --git a/bgzf.h b/bgzf.h
deleted file mode 100644
index 7295f37..0000000
--- a/bgzf.h
+++ /dev/null
@@ -1,157 +0,0 @@
-/* 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
-
-#include <stdint.h>
-#include <stdio.h>
-#include <zlib.h>
-#ifdef _USE_KNETFILE
-#include "knetfile.h"
-#endif
-
-//typedef int8_t bool;
-
-typedef struct {
- int file_descriptor;
- char open_mode; // 'r' or 'w'
- int16_t owned_file, compress_level;
-#ifdef _USE_KNETFILE
- union {
- knetFile *fpr;
- FILE *fpw;
- } x;
-#else
- FILE* file;
-#endif
- int uncompressed_block_size;
- int compressed_block_size;
- void* uncompressed_block;
- void* compressed_block;
- int64_t block_address;
- int block_length;
- int block_offset;
- int cache_size;
- const char* error;
- void *cache; // a pointer to a hash table
-} BGZF;
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-/*
- * Open an existing file descriptor for reading or writing.
- * Mode must be either "r" or "w".
- * A subsequent bgzf_close will not close the file descriptor.
- * Returns null on error.
- */
-BGZF* bgzf_fdopen(int fd, const char* __restrict mode);
-
-/*
- * Open the specified file for reading or writing.
- * Mode must be either "r" or "w".
- * Returns null on error.
- */
-BGZF* bgzf_open(const char* path, const char* __restrict mode);
-
-/*
- * Close the BGZ file and free all associated resources.
- * Does not close the underlying file descriptor if created with bgzf_fdopen.
- * Returns zero on success, -1 on error.
- */
-int bgzf_close(BGZF* fp);
-
-/*
- * Read up to length bytes from the file storing into data.
- * Returns the number of bytes actually read.
- * Returns zero on end of file.
- * Returns -1 on error.
- */
-int bgzf_read(BGZF* fp, void* data, int length);
-
-/*
- * Write length bytes from data to the file.
- * Returns the number of bytes written.
- * Returns -1 on error.
- */
-int bgzf_write(BGZF* fp, const void* data, int length);
-
-/*
- * Return a virtual file pointer to the current location in the file.
- * No interpetation of the value should be made, other than a subsequent
- * call to bgzf_seek can be used to position the file at the same point.
- * Return value is non-negative on success.
- * Returns -1 on error.
- */
-#define bgzf_tell(fp) ((fp->block_address << 16) | (fp->block_offset & 0xFFFF))
-
-/*
- * Set the file to read from the location specified by pos, which must
- * be a value previously returned by bgzf_tell for this file (but not
- * necessarily one returned by this file handle).
- * The where argument must be SEEK_SET.
- * Seeking on a file opened for write is not supported.
- * Returns zero on success, -1 on error.
- */
-int64_t bgzf_seek(BGZF* fp, int64_t pos, int where);
-
-/*
- * Set the cache size. Zero to disable. By default, caching is
- * disabled. The recommended cache size for frequent random access is
- * about 8M bytes.
- */
-void bgzf_set_cache_size(BGZF *fp, int cache_size);
-
-int bgzf_check_EOF(BGZF *fp);
-int bgzf_read_block(BGZF* fp);
-int bgzf_flush(BGZF* fp);
-int bgzf_flush_try(BGZF *fp, int size);
-int bgzf_check_bgzf(const char *fn);
-
-#ifdef __cplusplus
-}
-#endif
-
-static inline int bgzf_getc(BGZF *fp)
-{
- int c;
- if (fp->block_offset >= fp->block_length) {
- if (bgzf_read_block(fp) != 0) return -2; /* error */
- if (fp->block_length == 0) return -1; /* end-of-file */
- }
- c = ((unsigned char*)fp->uncompressed_block)[fp->block_offset++];
- if (fp->block_offset == fp->block_length) {
-#ifdef _USE_KNETFILE
- fp->block_address = knet_tell(fp->x.fpr);
-#else
- fp->block_address = ftello(fp->file);
-#endif
- fp->block_offset = 0;
- fp->block_length = 0;
- }
- return c;
-}
-
-#endif
diff --git a/bgzip.c b/bgzip.c
deleted file mode 100644
index ebcafa2..0000000
--- a/bgzip.c
+++ /dev/null
@@ -1,206 +0,0 @@
-/* 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 <fcntl.h>
-#include <unistd.h>
-#include <errno.h>
-#include <sys/select.h>
-#include <sys/stat.h>
-#include "bgzf.h"
-
-static const int WINDOW_SIZE = 64 * 1024;
-
-static int bgzip_main_usage()
-{
- fprintf(stderr, "\n");
- fprintf(stderr, "Usage: bgzip [options] [file] ...\n\n");
- fprintf(stderr, "Options: -c write on standard output, keep original files unchanged\n");
- fprintf(stderr, " -d decompress\n");
- fprintf(stderr, " -f overwrite files without asking\n");
- fprintf(stderr, " -b INT decompress at virtual file pointer INT\n");
- fprintf(stderr, " -s INT decompress INT bytes in the uncompressed file\n");
- fprintf(stderr, " -h give this help\n");
- fprintf(stderr, "\n");
- return 1;
-}
-
-static int write_open(const char *fn, int is_forced)
-{
- int fd = -1;
- char c;
- if (!is_forced) {
- if ((fd = open(fn, O_WRONLY | O_CREAT | O_TRUNC | O_EXCL, 0666)) < 0 && errno == EEXIST) {
- fprintf(stderr, "[bgzip] %s already exists; do you wish to overwrite (y or n)? ", fn);
- scanf("%c", &c);
- if (c != 'Y' && c != 'y') {
- fprintf(stderr, "[bgzip] not overwritten\n");
- exit(1);
- }
- }
- }
- if (fd < 0) {
- if ((fd = open(fn, O_WRONLY | O_CREAT | O_TRUNC, 0666)) < 0) {
- fprintf(stderr, "[bgzip] %s: Fail to write\n", fn);
- exit(1);
- }
- }
- return fd;
-}
-
-static void fail(BGZF* fp)
-{
- fprintf(stderr, "Error: %s\n", fp->error);
- exit(1);
-}
-
-int main(int argc, char **argv)
-{
- int c, compress, pstdout, is_forced;
- BGZF *fp;
- void *buffer;
- long start, end, size;
-
- compress = 1; pstdout = 0; start = 0; size = -1; end = -1; is_forced = 0;
- while((c = getopt(argc, argv, "cdhfb:s:")) >= 0){
- switch(c){
- case 'h': return bgzip_main_usage();
- case 'd': compress = 0; break;
- case 'c': pstdout = 1; break;
- case 'b': start = atol(optarg); break;
- case 's': size = atol(optarg); break;
- case 'f': is_forced = 1; break;
- }
- }
- if (size >= 0) end = start + size;
- if (end >= 0 && end < start) {
- fprintf(stderr, "[bgzip] Illegal region: [%ld, %ld]\n", start, end);
- return 1;
- }
- if (compress == 1) {
- struct stat sbuf;
- int f_src = fileno(stdin);
- int f_dst = fileno(stdout);
-
- if ( argc>optind )
- {
- if ( stat(argv[optind],&sbuf)<0 )
- {
- fprintf(stderr, "[bgzip] %s: %s\n", strerror(errno), argv[optind]);
- return 1;
- }
-
- if ((f_src = open(argv[optind], O_RDONLY)) < 0) {
- fprintf(stderr, "[bgzip] %s: %s\n", strerror(errno), argv[optind]);
- return 1;
- }
-
- if (pstdout)
- f_dst = fileno(stdout);
- else
- {
- char *name = malloc(strlen(argv[optind]) + 5);
- strcpy(name, argv[optind]);
- strcat(name, ".gz");
- f_dst = write_open(name, is_forced);
- if (f_dst < 0) return 1;
- free(name);
- }
- }
- else if (!pstdout && isatty(fileno((FILE *)stdout)) )
- return bgzip_main_usage();
-
- fp = bgzf_fdopen(f_dst, "w");
- buffer = malloc(WINDOW_SIZE);
- while ((c = read(f_src, buffer, WINDOW_SIZE)) > 0)
- if (bgzf_write(fp, buffer, c) < 0) fail(fp);
- // f_dst will be closed here
- if (bgzf_close(fp) < 0) fail(fp);
- if (argc > optind && !pstdout) unlink(argv[optind]);
- free(buffer);
- close(f_src);
- return 0;
- } else {
- struct stat sbuf;
- int f_dst;
-
- if ( argc>optind )
- {
- if ( stat(argv[optind],&sbuf)<0 )
- {
- fprintf(stderr, "[bgzip] %s: %s\n", strerror(errno), argv[optind]);
- return 1;
- }
- char *name;
- int len = strlen(argv[optind]);
- if ( strcmp(argv[optind]+len-3,".gz") )
- {
- fprintf(stderr, "[bgzip] %s: unknown suffix -- ignored\n", argv[optind]);
- return 1;
- }
- fp = bgzf_open(argv[optind], "r");
- if (fp == NULL) {
- fprintf(stderr, "[bgzip] Could not open file: %s\n", argv[optind]);
- return 1;
- }
-
- if (pstdout) {
- f_dst = fileno(stdout);
- }
- else {
- name = strdup(argv[optind]);
- name[strlen(name) - 3] = '\0';
- f_dst = write_open(name, is_forced);
- free(name);
- }
- }
- else if (!pstdout && isatty(fileno((FILE *)stdin)) )
- return bgzip_main_usage();
- else
- {
- f_dst = fileno(stdout);
- fp = bgzf_fdopen(fileno(stdin), "r");
- if (fp == NULL) {
- fprintf(stderr, "[bgzip] Could not read from stdin: %s\n", strerror(errno));
- return 1;
- }
- }
- buffer = malloc(WINDOW_SIZE);
- if (bgzf_seek(fp, start, SEEK_SET) < 0) fail(fp);
- while (1) {
- if (end < 0) c = bgzf_read(fp, buffer, WINDOW_SIZE);
- else c = bgzf_read(fp, buffer, (end - start > WINDOW_SIZE)? WINDOW_SIZE:(end - start));
- if (c == 0) break;
- if (c < 0) fail(fp);
- start += c;
- write(f_dst, buffer, c);
- if (end >= 0 && start >= end) break;
- }
- free(buffer);
- if (bgzf_close(fp) < 0) fail(fp);
- if (!pstdout) unlink(argv[optind]);
- return 0;
- }
-}
diff --git a/example.gtf.gz b/example.gtf.gz
deleted file mode 100644
index 693db0c..0000000
Binary files a/example.gtf.gz and /dev/null differ
diff --git a/example.gtf.gz.tbi b/example.gtf.gz.tbi
deleted file mode 100644
index 84a3c84..0000000
Binary files a/example.gtf.gz.tbi and /dev/null differ
diff --git a/index.c b/index.c
deleted file mode 100644
index 8f254e3..0000000
--- a/index.c
+++ /dev/null
@@ -1,998 +0,0 @@
-#include <ctype.h>
-#include <assert.h>
-#include <sys/stat.h>
-#include "khash.h"
-#include "ksort.h"
-#include "kstring.h"
-#include "bam_endian.h"
-#ifdef _USE_KNETFILE
-#include "knetfile.h"
-#endif
-#include "tabix.h"
-
-#define TAD_MIN_CHUNK_GAP 32768
-// 1<<14 is the size of minimum bin.
-#define TAD_LIDX_SHIFT 14
-
-typedef struct {
- uint64_t u, v;
-} pair64_t;
-
-#define pair64_lt(a,b) ((a).u < (b).u)
-KSORT_INIT(offt, pair64_t, pair64_lt)
-
-typedef struct {
- uint32_t m, n;
- pair64_t *list;
-} ti_binlist_t;
-
-typedef struct {
- int32_t n, m;
- uint64_t *offset;
-} ti_lidx_t;
-
-KHASH_MAP_INIT_INT(i, ti_binlist_t)
-KHASH_MAP_INIT_STR(s, int)
-
-struct __ti_index_t {
- ti_conf_t conf;
- int32_t n, max;
- khash_t(s) *tname;
- khash_t(i) **index;
- ti_lidx_t *index2;
-};
-
-struct __ti_iter_t {
- int from_first; // read from the first record; no random access
- int tid, beg, end, n_off, i, finished;
- uint64_t curr_off;
- kstring_t str;
- const ti_index_t *idx;
- pair64_t *off;
-};
-
-typedef struct {
- int tid, beg, end, bin;
-} ti_intv_t;
-
-ti_conf_t ti_conf_gff = { 0, 1, 4, 5, '#', 0 };
-ti_conf_t ti_conf_bed = { TI_FLAG_UCSC, 1, 2, 3, '#', 0 };
-ti_conf_t ti_conf_psltbl = { TI_FLAG_UCSC, 15, 17, 18, '#', 0 };
-ti_conf_t ti_conf_sam = { TI_PRESET_SAM, 3, 4, 0, '@', 0 };
-ti_conf_t ti_conf_vcf = { TI_PRESET_VCF, 1, 2, 0, '#', 0 };
-
-/***************
- * read a line *
- ***************/
-
-/*
-int ti_readline(BGZF *fp, kstring_t *str)
-{
- int c, l = 0;
- str->l = 0;
- while ((c = bgzf_getc(fp)) >= 0 && c != '\n') {
- ++l;
- if (c != '\r') kputc(c, str);
- }
- if (c < 0 && l == 0) return -1; // end of file
- return str->l;
-}
-*/
-
-/* Below is a faster implementation largely equivalent to the one
- * commented out above. */
-int ti_readline(BGZF *fp, kstring_t *str)
-{
- int l, state = 0;
- unsigned char *buf = (unsigned char*)fp->uncompressed_block;
- str->l = 0;
- do {
- if (fp->block_offset >= fp->block_length) {
- if (bgzf_read_block(fp) != 0) { state = -2; break; }
- if (fp->block_length == 0) { state = -1; break; }
- }
- for (l = fp->block_offset; l < fp->block_length && buf[l] != '\n'; ++l);
- if (l < fp->block_length) state = 1;
- l -= fp->block_offset;
- if (str->l + l + 1 >= str->m) {
- str->m = str->l + l + 2;
- kroundup32(str->m);
- str->s = (char*)realloc(str->s, str->m);
- }
- memcpy(str->s + str->l, buf + fp->block_offset, l);
- str->l += l;
- fp->block_offset += l + 1;
- if (fp->block_offset >= fp->block_length) {
-#ifdef _USE_KNETFILE
- fp->block_address = knet_tell(fp->x.fpr);
-#else
- fp->block_address = ftello(fp->file);
-#endif
- fp->block_offset = 0;
- fp->block_length = 0;
- }
- } while (state == 0);
- if (str->l == 0 && state < 0) return state;
- str->s[str->l] = 0;
- return str->l;
-}
-
-/*************************************
- * get the interval from a data line *
- *************************************/
-
-static inline int ti_reg2bin(uint32_t beg, uint32_t end)
-{
- --end;
- if (beg>>14 == end>>14) return 4681 + (beg>>14);
- if (beg>>17 == end>>17) return 585 + (beg>>17);
- if (beg>>20 == end>>20) return 73 + (beg>>20);
- if (beg>>23 == end>>23) return 9 + (beg>>23);
- if (beg>>26 == end>>26) return 1 + (beg>>26);
- return 0;
-}
-
-static int get_tid(ti_index_t *idx, const char *ss)
-{
- khint_t k;
- int tid;
- k = kh_get(s, idx->tname, ss);
- if (k == kh_end(idx->tname)) { // a new target sequence
- int ret, size;
- // update idx->n, ->max, ->index and ->index2
- if (idx->n == idx->max) {
- idx->max = idx->max? idx->max<<1 : 8;
- idx->index = realloc(idx->index, idx->max * sizeof(void*));
- idx->index2 = realloc(idx->index2, idx->max * sizeof(ti_lidx_t));
- }
- memset(&idx->index2[idx->n], 0, sizeof(ti_lidx_t));
- idx->index[idx->n++] = kh_init(i);
- // update ->tname
- tid = size = kh_size(idx->tname);
- k = kh_put(s, idx->tname, strdup(ss), &ret);
- kh_value(idx->tname, k) = size;
- assert(idx->n == kh_size(idx->tname));
- } else tid = kh_value(idx->tname, k);
- return tid;
-}
-
-int ti_get_intv(const ti_conf_t *conf, int len, char *line, ti_interval_t *intv)
-{
- int i, b = 0, id = 1, ncols = 0;
- char *s;
- intv->ss = intv->se = 0; intv->beg = intv->end = -1;
- for (i = 0; i <= len; ++i) {
- if (line[i] == '\t' || line[i] == 0) {
- ++ncols;
- if (id == conf->sc) {
- intv->ss = line + b; intv->se = line + i;
- } else if (id == conf->bc) {
- // here ->beg is 0-based.
- intv->beg = intv->end = strtol(line + b, &s, 0);
- if (!(conf->preset&TI_FLAG_UCSC)) --intv->beg;
- else ++intv->end;
- if (intv->beg < 0) intv->beg = 0;
- if (intv->end < 1) intv->end = 1;
- } else {
- if ((conf->preset&0xffff) == TI_PRESET_GENERIC) {
- if (id == conf->ec) intv->end = strtol(line + b, &s, 0);
- } else if ((conf->preset&0xffff) == TI_PRESET_SAM) {
- if (id == 6) { // CIGAR
- int l = 0, op;
- char *t;
- for (s = line + b; s < line + i;) {
- long x = strtol(s, &t, 10);
- op = toupper(*t);
- if (op == 'M' || op == 'D' || op == 'N') l += x;
- s = t + 1;
- }
- if (l == 0) l = 1;
- intv->end = intv->beg + l;
- }
- } else if ((conf->preset&0xffff) == TI_PRESET_VCF) {
- // FIXME: the following is NOT tested and is likely to be buggy
- if (id == 4) {
- if (b < i) intv->end = intv->beg + (i - b);
- } else if (id == 8) { // look for "END="
- int c = line[i];
- line[i] = 0;
- s = strstr(line + b, "END=");
- if (s == line + b) s += 4;
- else if (s) {
- s = strstr(line + b, ";END=");
- if (s) s += 5;
- }
- if (s) intv->end = strtol(s, &s, 0);
- line[i] = c;
- }
- }
- }
- b = i + 1;
- ++id;
- }
- }
-/*
- if (ncols < conf->sc || ncols < conf->bc || ncols < conf->ec) {
- if (ncols == 1) fprintf(stderr,"[get_intv] Is the file tab-delimited? The line has %d field only: %s\n", ncols, line);
- else fprintf(stderr,"[get_intv] The line has %d field(s) only: %s\n", ncols, line);
- exit(1);
- }
-*/
- if (intv->ss == 0 || intv->se == 0 || intv->beg < 0 || intv->end < 0) return -1;
- return 0;
-}
-
-static int get_intv(ti_index_t *idx, kstring_t *str, ti_intv_t *intv)
-{
- ti_interval_t x;
- intv->tid = intv->beg = intv->end = intv->bin = -1;
- if (ti_get_intv(&idx->conf, str->l, str->s, &x) == 0) {
- int c = *x.se;
- *x.se = '\0'; intv->tid = get_tid(idx, x.ss); *x.se = c;
- intv->beg = x.beg; intv->end = x.end;
- intv->bin = ti_reg2bin(intv->beg, intv->end);
- return (intv->tid >= 0 && intv->beg >= 0 && intv->end >= 0)? 0 : -1;
- } else {
- fprintf(stderr, "[%s] the following line cannot be parsed and skipped: %s\n", __func__, str->s);
- return -1;
- }
-}
-
-/************
- * indexing *
- ************/
-
-// requirement: len <= LEN_MASK
-static inline void insert_offset(khash_t(i) *h, int bin, uint64_t beg, uint64_t end)
-{
- khint_t k;
- ti_binlist_t *l;
- int ret;
- k = kh_put(i, h, bin, &ret);
- l = &kh_value(h, k);
- if (ret) { // not present
- l->m = 1; l->n = 0;
- l->list = (pair64_t*)calloc(l->m, 16);
- }
- if (l->n == l->m) {
- l->m <<= 1;
- l->list = (pair64_t*)realloc(l->list, l->m * 16);
- }
- l->list[l->n].u = beg; l->list[l->n++].v = end;
-}
-
-static inline uint64_t insert_offset2(ti_lidx_t *index2, int _beg, int _end, uint64_t offset)
-{
- int i, beg, end;
- beg = _beg >> TAD_LIDX_SHIFT;
- end = (_end - 1) >> TAD_LIDX_SHIFT;
- if (index2->m < end + 1) {
- int old_m = index2->m;
- index2->m = end + 1;
- kroundup32(index2->m);
- index2->offset = (uint64_t*)realloc(index2->offset, index2->m * 8);
- memset(index2->offset + old_m, 0, 8 * (index2->m - old_m));
- }
- if (beg == end) {
- if (index2->offset[beg] == 0) index2->offset[beg] = offset;
- } else {
- for (i = beg; i <= end; ++i)
- if (index2->offset[i] == 0) index2->offset[i] = offset;
- }
- if (index2->n < end + 1) index2->n = end + 1;
- return (uint64_t)beg<<32 | end;
-}
-
-static void merge_chunks(ti_index_t *idx)
-{
- khash_t(i) *index;
- int i, l, m;
- khint_t k;
- for (i = 0; i < idx->n; ++i) {
- index = idx->index[i];
- for (k = kh_begin(index); k != kh_end(index); ++k) {
- ti_binlist_t *p;
- if (!kh_exist(index, k)) continue;
- p = &kh_value(index, k);
- m = 0;
- for (l = 1; l < p->n; ++l) {
- if (p->list[m].v>>16 == p->list[l].u>>16) p->list[m].v = p->list[l].v;
- else p->list[++m] = p->list[l];
- } // ~for(l)
- p->n = m + 1;
- } // ~for(k)
- } // ~for(i)
-}
-
-static void fill_missing(ti_index_t *idx)
-{
- int i, j;
- for (i = 0; i < idx->n; ++i) {
- ti_lidx_t *idx2 = &idx->index2[i];
- for (j = 1; j < idx2->n; ++j)
- if (idx2->offset[j] == 0)
- idx2->offset[j] = idx2->offset[j-1];
- }
-}
-
-ti_index_t *ti_index_core(BGZF *fp, const ti_conf_t *conf)
-{
- int ret;
- ti_index_t *idx;
- uint32_t last_bin, save_bin;
- int32_t last_coor, last_tid, save_tid;
- uint64_t save_off, last_off, lineno = 0, offset0 = (uint64_t)-1, tmp;
- kstring_t *str;
-
- str = calloc(1, sizeof(kstring_t));
-
- idx = (ti_index_t*)calloc(1, sizeof(ti_index_t));
- idx->conf = *conf;
- idx->n = idx->max = 0;
- idx->tname = kh_init(s);
- idx->index = 0;
- idx->index2 = 0;
-
- save_bin = save_tid = last_tid = last_bin = 0xffffffffu;
- save_off = last_off = bgzf_tell(fp); last_coor = 0xffffffffu;
- while ((ret = ti_readline(fp, str)) >= 0) {
- ti_intv_t intv;
- ++lineno;
- if (lineno <= idx->conf.line_skip || str->s[0] == idx->conf.meta_char) {
- last_off = bgzf_tell(fp);
- continue;
- }
- get_intv(idx, str, &intv);
- if ( intv.beg<0 || intv.end<0 )
- {
- fprintf(stderr,"[ti_index_core] the indexes overlap or are out of bounds\n");
- exit(1);
- }
- if (last_tid != intv.tid) { // change of chromosomes
- if (last_tid>intv.tid )
- {
- fprintf(stderr,"[ti_index_core] the chromosome blocks not continuous at line %llu, is the file sorted? [pos %d]\n",(unsigned long long)lineno,intv.beg+1);
- exit(1);
- }
- last_tid = intv.tid;
- last_bin = 0xffffffffu;
- } else if (last_coor > intv.beg) {
- fprintf(stderr, "[ti_index_core] the file out of order at line %llu\n", (unsigned long long)lineno);
- exit(1);
- }
- tmp = insert_offset2(&idx->index2[intv.tid], intv.beg, intv.end, last_off);
- if (last_off == 0) offset0 = tmp;
- if (intv.bin != last_bin) { // then possibly write the binning index
- if (save_bin != 0xffffffffu) // save_bin==0xffffffffu only happens to the first record
- insert_offset(idx->index[save_tid], save_bin, save_off, last_off);
- save_off = last_off;
- save_bin = last_bin = intv.bin;
- save_tid = intv.tid;
- if (save_tid < 0) break;
- }
- if (bgzf_tell(fp) <= last_off) {
- fprintf(stderr, "[ti_index_core] bug in BGZF: %llx < %llx\n",
- (unsigned long long)bgzf_tell(fp), (unsigned long long)last_off);
- exit(1);
- }
- last_off = bgzf_tell(fp);
- last_coor = intv.beg;
- }
- if (save_tid >= 0) insert_offset(idx->index[save_tid], save_bin, save_off, bgzf_tell(fp));
- merge_chunks(idx);
- fill_missing(idx);
- if (offset0 != (uint64_t)-1 && idx->n && idx->index2[0].offset) {
- int i, beg = offset0>>32, end = offset0&0xffffffffu;
- for (i = beg; i <= end; ++i) idx->index2[0].offset[i] = 0;
- }
-
- free(str->s); free(str);
- return idx;
-}
-
-void ti_index_destroy(ti_index_t *idx)
-{
- khint_t k;
- int i;
- if (idx == 0) return;
- // destroy the name hash table
- for (k = kh_begin(idx->tname); k != kh_end(idx->tname); ++k) {
- if (kh_exist(idx->tname, k))
- free((char*)kh_key(idx->tname, k));
- }
- kh_destroy(s, idx->tname);
- // destroy the binning index
- for (i = 0; i < idx->n; ++i) {
- khash_t(i) *index = idx->index[i];
- ti_lidx_t *index2 = idx->index2 + i;
- for (k = kh_begin(index); k != kh_end(index); ++k) {
- if (kh_exist(index, k))
- free(kh_value(index, k).list);
- }
- kh_destroy(i, index);
- free(index2->offset);
- }
- free(idx->index);
- // destroy the linear index
- free(idx->index2);
- free(idx);
-}
-
-/******************
- * index file I/O *
- ******************/
-
-void ti_index_save(const ti_index_t *idx, BGZF *fp)
-{
- int32_t i, size, ti_is_be;
- khint_t k;
- ti_is_be = bam_is_big_endian();
- bgzf_write(fp, "TBI\1", 4);
- if (ti_is_be) {
- uint32_t x = idx->n;
- bgzf_write(fp, bam_swap_endian_4p(&x), 4);
- } else bgzf_write(fp, &idx->n, 4);
- assert(sizeof(ti_conf_t) == 24);
- if (ti_is_be) { // write ti_conf_t;
- uint32_t x[6];
- memcpy(x, &idx->conf, 24);
- for (i = 0; i < 6; ++i) bgzf_write(fp, bam_swap_endian_4p(&x[i]), 4);
- } else bgzf_write(fp, &idx->conf, sizeof(ti_conf_t));
- { // write target names
- char **name;
- int32_t l = 0;
- name = calloc(kh_size(idx->tname), sizeof(void*));
- for (k = kh_begin(idx->tname); k != kh_end(idx->tname); ++k)
- if (kh_exist(idx->tname, k))
- name[kh_value(idx->tname, k)] = (char*)kh_key(idx->tname, k);
- for (i = 0; i < kh_size(idx->tname); ++i)
- l += strlen(name[i]) + 1;
- if (ti_is_be) bgzf_write(fp, bam_swap_endian_4p(&l), 4);
- else bgzf_write(fp, &l, 4);
- for (i = 0; i < kh_size(idx->tname); ++i)
- bgzf_write(fp, name[i], strlen(name[i]) + 1);
- free(name);
- }
- for (i = 0; i < idx->n; ++i) {
- khash_t(i) *index = idx->index[i];
- ti_lidx_t *index2 = idx->index2 + i;
- // write binning index
- size = kh_size(index);
- if (ti_is_be) { // big endian
- uint32_t x = size;
- bgzf_write(fp, bam_swap_endian_4p(&x), 4);
- } else bgzf_write(fp, &size, 4);
- for (k = kh_begin(index); k != kh_end(index); ++k) {
- if (kh_exist(index, k)) {
- ti_binlist_t *p = &kh_value(index, k);
- if (ti_is_be) { // big endian
- uint32_t x;
- x = kh_key(index, k); bgzf_write(fp, bam_swap_endian_4p(&x), 4);
- x = p->n; bgzf_write(fp, bam_swap_endian_4p(&x), 4);
- for (x = 0; (int)x < p->n; ++x) {
- bam_swap_endian_8p(&p->list[x].u);
- bam_swap_endian_8p(&p->list[x].v);
- }
- bgzf_write(fp, p->list, 16 * p->n);
- for (x = 0; (int)x < p->n; ++x) {
- bam_swap_endian_8p(&p->list[x].u);
- bam_swap_endian_8p(&p->list[x].v);
- }
- } else {
- bgzf_write(fp, &kh_key(index, k), 4);
- bgzf_write(fp, &p->n, 4);
- bgzf_write(fp, p->list, 16 * p->n);
- }
- }
- }
- // write linear index (index2)
- if (ti_is_be) {
- int x = index2->n;
- bgzf_write(fp, bam_swap_endian_4p(&x), 4);
- } else bgzf_write(fp, &index2->n, 4);
- if (ti_is_be) { // big endian
- int x;
- for (x = 0; (int)x < index2->n; ++x)
- bam_swap_endian_8p(&index2->offset[x]);
- bgzf_write(fp, index2->offset, 8 * index2->n);
- for (x = 0; (int)x < index2->n; ++x)
- bam_swap_endian_8p(&index2->offset[x]);
- } else bgzf_write(fp, index2->offset, 8 * index2->n);
- }
-}
-
-static ti_index_t *ti_index_load_core(BGZF *fp)
-{
- int i, ti_is_be;
- char magic[4];
- ti_index_t *idx;
- ti_is_be = bam_is_big_endian();
- if (fp == 0) {
- fprintf(stderr, "[ti_index_load_core] fail to load index.\n");
- return 0;
- }
- bgzf_read(fp, magic, 4);
- if (strncmp(magic, "TBI\1", 4)) {
- fprintf(stderr, "[ti_index_load] wrong magic number.\n");
- return 0;
- }
- idx = (ti_index_t*)calloc(1, sizeof(ti_index_t));
- bgzf_read(fp, &idx->n, 4);
- if (ti_is_be) bam_swap_endian_4p(&idx->n);
- idx->tname = kh_init(s);
- idx->index = (khash_t(i)**)calloc(idx->n, sizeof(void*));
- idx->index2 = (ti_lidx_t*)calloc(idx->n, sizeof(ti_lidx_t));
- // read idx->conf
- bgzf_read(fp, &idx->conf, sizeof(ti_conf_t));
- if (ti_is_be) {
- bam_swap_endian_4p(&idx->conf.preset);
- bam_swap_endian_4p(&idx->conf.sc);
- bam_swap_endian_4p(&idx->conf.bc);
- bam_swap_endian_4p(&idx->conf.ec);
- bam_swap_endian_4p(&idx->conf.meta_char);
- bam_swap_endian_4p(&idx->conf.line_skip);
- }
- { // read target names
- int j, ret;
- kstring_t *str;
- int32_t l;
- uint8_t *buf;
- bgzf_read(fp, &l, 4);
- if (ti_is_be) bam_swap_endian_4p(&l);
- buf = calloc(l, 1);
- bgzf_read(fp, buf, l);
- str = calloc(1, sizeof(kstring_t));
- for (i = j = 0; i < l; ++i) {
- if (buf[i] == 0) {
- khint_t k = kh_put(s, idx->tname, strdup(str->s), &ret);
- kh_value(idx->tname, k) = j++;
- str->l = 0;
- } else kputc(buf[i], str);
- }
- free(str->s); free(str); free(buf);
- }
- for (i = 0; i < idx->n; ++i) {
- khash_t(i) *index;
- ti_lidx_t *index2 = idx->index2 + i;
- uint32_t key, size;
- khint_t k;
- int j, ret;
- ti_binlist_t *p;
- index = idx->index[i] = kh_init(i);
- // load binning index
- bgzf_read(fp, &size, 4);
- if (ti_is_be) bam_swap_endian_4p(&size);
- for (j = 0; j < (int)size; ++j) {
- bgzf_read(fp, &key, 4);
- if (ti_is_be) bam_swap_endian_4p(&key);
- k = kh_put(i, index, key, &ret);
- p = &kh_value(index, k);
- bgzf_read(fp, &p->n, 4);
- if (ti_is_be) bam_swap_endian_4p(&p->n);
- p->m = p->n;
- p->list = (pair64_t*)malloc(p->m * 16);
- bgzf_read(fp, p->list, 16 * p->n);
- if (ti_is_be) {
- int x;
- for (x = 0; x < p->n; ++x) {
- bam_swap_endian_8p(&p->list[x].u);
- bam_swap_endian_8p(&p->list[x].v);
- }
- }
- }
- // load linear index
- bgzf_read(fp, &index2->n, 4);
- if (ti_is_be) bam_swap_endian_4p(&index2->n);
- index2->m = index2->n;
- index2->offset = (uint64_t*)calloc(index2->m, 8);
- bgzf_read(fp, index2->offset, index2->n * 8);
- if (ti_is_be)
- for (j = 0; j < index2->n; ++j) bam_swap_endian_8p(&index2->offset[j]);
- }
- return idx;
-}
-
-ti_index_t *ti_index_load_local(const char *fnidx)
-{
- BGZF *fp;
- fp = bgzf_open(fnidx, "r");
- if (fp) {
- ti_index_t *idx = ti_index_load_core(fp);
- bgzf_close(fp);
- return idx;
- } else return 0;
-}
-
-#ifdef _USE_KNETFILE
-static void download_from_remote(const char *url)
-{
- const int buf_size = 1 * 1024 * 1024;
- char *fn;
- FILE *fp;
- uint8_t *buf;
- knetFile *fp_remote;
- int l;
- 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;
- ++fn; // fn now points to the file name
- fp_remote = knet_open(url, "r");
- if (fp_remote == 0) {
- fprintf(stderr, "[download_from_remote] fail to open remote file.\n");
- return;
- }
- if ((fp = fopen(fn, "w")) == 0) {
- fprintf(stderr, "[download_from_remote] fail to create file in the working directory.\n");
- knet_close(fp_remote);
- return;
- }
- buf = (uint8_t*)calloc(buf_size, 1);
- while ((l = knet_read(fp_remote, buf, buf_size)) != 0)
- fwrite(buf, 1, l, fp);
- free(buf);
- fclose(fp);
- knet_close(fp_remote);
-}
-#else
-static void download_from_remote(const char *url)
-{
- return;
-}
-#endif
-
-static char *get_local_version(const char *fn)
-{
- struct stat sbuf;
- char *fnidx = (char*)calloc(strlen(fn) + 5, 1);
- strcat(strcpy(fnidx, fn), ".tbi");
- if ((strstr(fnidx, "ftp://") == fnidx || strstr(fnidx, "http://") == fnidx)) {
- char *p, *url;
- int l = strlen(fnidx);
- for (p = fnidx + l - 1; p >= fnidx; --p)
- if (*p == '/') break;
- url = fnidx; fnidx = strdup(p + 1);
- if (stat(fnidx, &sbuf) == 0) {
- free(url);
- return fnidx;
- }
- fprintf(stderr, "[%s] downloading the index file...\n", __func__);
- download_from_remote(url);
- free(url);
- }
- if (stat(fnidx, &sbuf) == 0) return fnidx;
- free(fnidx); return 0;
-}
-
-const char **ti_seqname(const ti_index_t *idx, int *n)
-{
- const char **names;
- khint_t k;
- *n = idx->n;
- names = calloc(idx->n, sizeof(void*));
- for (k = kh_begin(idx->tname); k < kh_end(idx->tname); ++k)
- if (kh_exist(idx->tname, k))
- names[kh_val(idx->tname, k)] = kh_key(idx->tname, k);
- return names;
-}
-
-ti_index_t *ti_index_load(const char *fn)
-{
- ti_index_t *idx;
- char *fname = get_local_version(fn);
- if (fname == 0) return 0;
- idx = ti_index_load_local(fname);
- if (idx == 0) fprintf(stderr, "[ti_index_load] fail to load the index: %s\n", fname);
- free(fname);
- return idx;
-}
-
-int ti_index_build2(const char *fn, const ti_conf_t *conf, const char *_fnidx)
-{
- char *fnidx;
- BGZF *fp, *fpidx;
- ti_index_t *idx;
- if ((fp = bgzf_open(fn, "r")) == 0) {
- fprintf(stderr, "[ti_index_build2] fail to open the file: %s\n", fn);
- return -1;
- }
- idx = ti_index_core(fp, conf);
- bgzf_close(fp);
- if (_fnidx == 0) {
- fnidx = (char*)calloc(strlen(fn) + 5, 1);
- strcpy(fnidx, fn); strcat(fnidx, ".tbi");
- } else fnidx = strdup(_fnidx);
- fpidx = bgzf_open(fnidx, "w");
- if (fpidx == 0) {
- fprintf(stderr, "[ti_index_build2] fail to create the index file.\n");
- free(fnidx);
- return -1;
- }
- ti_index_save(idx, fpidx);
- ti_index_destroy(idx);
- bgzf_close(fpidx);
- free(fnidx);
- return 0;
-}
-
-int ti_index_build(const char *fn, const ti_conf_t *conf)
-{
- return ti_index_build2(fn, conf, 0);
-}
-
-/********************************************
- * parse a region in the format chr:beg-end *
- ********************************************/
-
-int ti_get_tid(const ti_index_t *idx, const char *name)
-{
- khiter_t iter;
- const khash_t(s) *h = idx->tname;
- iter = kh_get(s, h, name); /* get the tid */
- if (iter == kh_end(h)) return -1;
- return kh_value(h, iter);
-}
-
-int ti_parse_region(const ti_index_t *idx, const char *str, int *tid, int *begin, int *end)
-{
- char *s, *p;
- int i, l, k;
- l = strlen(str);
- p = s = (char*)malloc(l+1);
- /* squeeze out "," */
- for (i = k = 0; i != l; ++i)
- if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
- s[k] = 0;
- for (i = 0; i != k; ++i) if (s[i] == ':') break;
- s[i] = 0;
- if ((*tid = ti_get_tid(idx, s)) < 0) {
- free(s);
- return -1;
- }
- if (i == k) { /* dump the whole sequence */
- *begin = 0; *end = 1<<29; free(s);
- return 0;
- }
- for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
- *begin = atoi(p);
- if (i < k) {
- p = s + i + 1;
- *end = atoi(p);
- } else *end = 1<<29;
- if (*begin > 0) --*begin;
- free(s);
- if (*begin > *end) return -1;
- return 0;
-}
-
-/*******************************
- * retrieve a specified region *
- *******************************/
-
-#define MAX_BIN 37450 // =(8^6-1)/7+1
-
-static inline int reg2bins(uint32_t beg, uint32_t end, uint16_t list[MAX_BIN])
-{
- int i = 0, k;
- if (beg >= end) return 0;
- if (end >= 1u<<29) end = 1u<<29;
- --end;
- list[i++] = 0;
- for (k = 1 + (beg>>26); k <= 1 + (end>>26); ++k) list[i++] = k;
- for (k = 9 + (beg>>23); k <= 9 + (end>>23); ++k) list[i++] = k;
- for (k = 73 + (beg>>20); k <= 73 + (end>>20); ++k) list[i++] = k;
- for (k = 585 + (beg>>17); k <= 585 + (end>>17); ++k) list[i++] = k;
- for (k = 4681 + (beg>>14); k <= 4681 + (end>>14); ++k) list[i++] = k;
- return i;
-}
-
-ti_iter_t ti_iter_first()
-{
- ti_iter_t iter;
- iter = calloc(1, sizeof(struct __ti_iter_t));
- iter->from_first = 1;
- return iter;
-}
-
-ti_iter_t ti_iter_query(const ti_index_t *idx, int tid, int beg, int end)
-{
- uint16_t *bins;
- int i, n_bins, n_off;
- pair64_t *off;
- khint_t k;
- khash_t(i) *index;
- uint64_t min_off;
- ti_iter_t iter = 0;
-
- if (beg < 0) beg = 0;
- if (end < beg) return 0;
- // initialize the iterator
- iter = calloc(1, sizeof(struct __ti_iter_t));
- iter->idx = idx; iter->tid = tid; iter->beg = beg; iter->end = end; iter->i = -1;
- // random access
- bins = (uint16_t*)calloc(MAX_BIN, 2);
- n_bins = reg2bins(beg, end, bins);
- index = idx->index[tid];
- if (idx->index2[tid].n > 0) {
- min_off = (beg>>TAD_LIDX_SHIFT >= idx->index2[tid].n)? idx->index2[tid].offset[idx->index2[tid].n-1]
- : idx->index2[tid].offset[beg>>TAD_LIDX_SHIFT];
- if (min_off == 0) { // improvement for index files built by tabix prior to 0.1.4
- int n = beg>>TAD_LIDX_SHIFT;
- if (n > idx->index2[tid].n) n = idx->index2[tid].n;
- for (i = n - 1; i >= 0; --i)
- if (idx->index2[tid].offset[i] != 0) break;
- if (i >= 0) min_off = idx->index2[tid].offset[i];
- }
- } else min_off = 0; // tabix 0.1.2 may produce such index files
- for (i = n_off = 0; i < n_bins; ++i) {
- if ((k = kh_get(i, index, bins[i])) != kh_end(index))
- n_off += kh_value(index, k).n;
- }
- if (n_off == 0) {
- free(bins); return iter;
- }
- off = (pair64_t*)calloc(n_off, 16);
- for (i = n_off = 0; i < n_bins; ++i) {
- if ((k = kh_get(i, index, bins[i])) != kh_end(index)) {
- int j;
- ti_binlist_t *p = &kh_value(index, k);
- for (j = 0; j < p->n; ++j)
- if (p->list[j].v > min_off) off[n_off++] = p->list[j];
- }
- }
- if (n_off == 0) {
- free(bins); free(off); return iter;
- }
- free(bins);
- {
- int l;
- ks_introsort(offt, n_off, off);
- // resolve completely contained adjacent blocks
- for (i = 1, l = 0; i < n_off; ++i)
- if (off[l].v < off[i].v)
- off[++l] = off[i];
- n_off = l + 1;
- // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
- for (i = 1; i < n_off; ++i)
- if (off[i-1].v >= off[i].u) off[i-1].v = off[i].u;
- { // merge adjacent blocks
- for (i = 1, l = 0; i < n_off; ++i) {
- if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v;
- else off[++l] = off[i];
- }
- n_off = l + 1;
- }
- }
- iter->n_off = n_off; iter->off = off;
- return iter;
-}
-
-const char *ti_iter_read(BGZF *fp, ti_iter_t iter, int *len)
-{
- if (iter->finished) return 0;
- if (iter->from_first) {
- int ret;
- if ((ret = ti_readline(fp, &iter->str)) < 0) {
- iter->finished = 1;
- return 0;
- } else {
- if (len) *len = iter->str.l;
- return iter->str.s;
- }
- }
- if (iter->n_off == 0) return 0;
- while (1) {
- int ret;
- if (iter->curr_off == 0 || iter->curr_off >= iter->off[iter->i].v) { // then jump to the next chunk
- if (iter->i == iter->n_off - 1) break; // no more chunks
- if (iter->i >= 0) assert(iter->curr_off == iter->off[iter->i].v); // otherwise bug
- if (iter->i < 0 || iter->off[iter->i].v != iter->off[iter->i+1].u) { // not adjacent chunks; then seek
- bgzf_seek(fp, iter->off[iter->i+1].u, SEEK_SET);
- iter->curr_off = bgzf_tell(fp);
- }
- ++iter->i;
- }
- if ((ret = ti_readline(fp, &iter->str)) >= 0) {
- ti_intv_t intv;
- iter->curr_off = bgzf_tell(fp);
- if (iter->str.s[0] == iter->idx->conf.meta_char) continue;
- get_intv((ti_index_t*)iter->idx, &iter->str, &intv);
- if (intv.tid != iter->tid || intv.beg >= iter->end) break; // no need to proceed
- else if (intv.end > iter->beg && iter->end > intv.beg) {
- if (len) *len = iter->str.l;
- return iter->str.s;
- }
- } else break; // end of file
- }
- iter->finished = 1;
- return 0;
-}
-
-void ti_iter_destroy(ti_iter_t iter)
-{
- if (iter) {
- free(iter->str.s); free(iter->off);
- free(iter);
- }
-}
-
-int ti_fetch(BGZF *fp, const ti_index_t *idx, int tid, int beg, int end, void *data, ti_fetch_f func)
-{
- ti_iter_t iter;
- const char *s;
- int len;
- iter = ti_iter_query(idx, tid, beg, end);
- while ((s = ti_iter_read(fp, iter, &len)) != 0)
- func(len, s, data);
- ti_iter_destroy(iter);
- return 0;
-}
-
-const ti_conf_t *ti_get_conf(ti_index_t *idx) { return idx? &idx->conf : 0; }
-
-/*******************
- * High-level APIs *
- *******************/
-
-tabix_t *ti_open(const char *fn, const char *fnidx)
-{
- tabix_t *t;
- BGZF *fp;
- if ((fp = bgzf_open(fn, "r")) == 0) return 0;
- t = calloc(1, sizeof(tabix_t));
- t->fn = strdup(fn);
- if (fnidx) t->fnidx = strdup(fnidx);
- t->fp = fp;
- return t;
-}
-
-void ti_close(tabix_t *t)
-{
- if (t) {
- bgzf_close(t->fp);
- if (t->idx) ti_index_destroy(t->idx);
- free(t->fn); free(t->fnidx);
- free(t);
- }
-}
-
-int ti_lazy_index_load(tabix_t *t)
-{
- if (t->idx == 0) { // load index
- if (t->fnidx) t->idx = ti_index_load_local(t->fnidx);
- else t->idx = ti_index_load(t->fn);
- if (t->idx == 0) return -1; // fail to load index
- }
- return 0;
-}
-
-ti_iter_t ti_queryi(tabix_t *t, int tid, int beg, int end)
-{
- if (tid < 0) return ti_iter_first();
- if (ti_lazy_index_load(t) != 0) return 0;
- return ti_iter_query(t->idx, tid, beg, end);
-}
-
-ti_iter_t ti_querys(tabix_t *t, const char *reg)
-{
- int tid, beg, end;
- if (reg == 0) return ti_iter_first();
- if (ti_lazy_index_load(t) != 0) return 0;
- if (ti_parse_region(t->idx, reg, &tid, &beg, &end) < 0) return 0;
- return ti_iter_query(t->idx, tid, beg, end);
-}
-
-ti_iter_t ti_query(tabix_t *t, const char *name, int beg, int end)
-{
- int tid;
- if (name == 0) return ti_iter_first();
- // then need to load the index
- if (ti_lazy_index_load(t) != 0) return 0;
- if ((tid = ti_get_tid(t->idx, name)) < 0) return 0;
- return ti_iter_query(t->idx, tid, beg, end);
-}
-
-const char *ti_read(tabix_t *t, ti_iter_t iter, int *len)
-{
- return ti_iter_read(t->fp, iter, len);
-}
diff --git a/khash.h b/khash.h
deleted file mode 100644
index 1d583ef..0000000
--- a/khash.h
+++ /dev/null
@@ -1,486 +0,0 @@
-/* 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 at sanger.ac.uk> */
-
-/*
- An example:
-
-#include "khash.h"
-KHASH_MAP_INIT_INT(32, char)
-int main() {
- int ret, is_missing;
- khiter_t k;
- khash_t(32) *h = kh_init(32);
- k = kh_put(32, h, 5, &ret);
- if (!ret) kh_del(32, h, k);
- kh_value(h, k) = 10;
- k = kh_get(32, h, 10);
- is_missing = (k == kh_end(h));
- k = kh_get(32, h, 5);
- kh_del(32, h, k);
- for (k = kh_begin(h); k != kh_end(h); ++k)
- if (kh_exist(h, k)) kh_value(h, k) = 1;
- kh_destroy(32, h);
- return 0;
-}
-*/
-
-/*
- 2008-09-19 (0.2.3):
-
- * Corrected the example
- * Improved interfaces
-
- 2008-09-11 (0.2.2):
-
- * Improved speed a little in kh_put()
-
- 2008-09-10 (0.2.1):
-
- * Added kh_clear()
- * Fixed a compiling error
-
- 2008-09-02 (0.2.0):
-
- * Changed to token concatenation which increases flexibility.
-
- 2008-08-31 (0.1.2):
-
- * Fixed a bug in kh_get(), which has not been tested previously.
-
- 2008-08-31 (0.1.1):
-
- * Added destructor
-*/
-
-
-#ifndef __AC_KHASH_H
-#define __AC_KHASH_H
-
-/*!
- @header
-
- Generic hash table library.
-
- @copyright Heng Li
- */
-
-#define AC_VERSION_KHASH_H "0.2.2"
-
-#include <stdint.h>
-#include <stdlib.h>
-#include <string.h>
-
-typedef uint32_t khint_t;
-typedef khint_t khiter_t;
-
-#define __ac_HASH_PRIME_SIZE 32
-static const uint32_t __ac_prime_list[__ac_HASH_PRIME_SIZE] =
-{
- 0ul, 3ul, 11ul, 23ul, 53ul,
- 97ul, 193ul, 389ul, 769ul, 1543ul,
- 3079ul, 6151ul, 12289ul, 24593ul, 49157ul,
- 98317ul, 196613ul, 393241ul, 786433ul, 1572869ul,
- 3145739ul, 6291469ul, 12582917ul, 25165843ul, 50331653ul,
- 100663319ul, 201326611ul, 402653189ul, 805306457ul, 1610612741ul,
- 3221225473ul, 4294967291ul
-};
-
-#define __ac_isempty(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&2)
-#define __ac_isdel(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&1)
-#define __ac_iseither(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&3)
-#define __ac_set_isdel_false(flag, i) (flag[i>>4]&=~(1ul<<((i&0xfU)<<1)))
-#define __ac_set_isempty_false(flag, i) (flag[i>>4]&=~(2ul<<((i&0xfU)<<1)))
-#define __ac_set_isboth_false(flag, i) (flag[i>>4]&=~(3ul<<((i&0xfU)<<1)))
-#define __ac_set_isdel_true(flag, i) (flag[i>>4]|=1ul<<((i&0xfU)<<1))
-
-static const double __ac_HASH_UPPER = 0.77;
-
-#define KHASH_INIT(name, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \
- typedef struct { \
- khint_t n_buckets, size, n_occupied, upper_bound; \
- uint32_t *flags; \
- khkey_t *keys; \
- khval_t *vals; \
- } kh_##name##_t; \
- static inline kh_##name##_t *kh_init_##name() { \
- return (kh_##name##_t*)calloc(1, sizeof(kh_##name##_t)); \
- } \
- static inline void kh_destroy_##name(kh_##name##_t *h) \
- { \
- if (h) { \
- free(h->keys); free(h->flags); \
- free(h->vals); \
- free(h); \
- } \
- } \
- static inline void kh_clear_##name(kh_##name##_t *h) \
- { \
- if (h && h->flags) { \
- memset(h->flags, 0xaa, ((h->n_buckets>>4) + 1) * sizeof(uint32_t)); \
- h->size = h->n_occupied = 0; \
- } \
- } \
- static inline khint_t kh_get_##name(const kh_##name##_t *h, khkey_t key) \
- { \
- if (h->n_buckets) { \
- khint_t inc, k, i, last; \
- k = __hash_func(key); i = k % h->n_buckets; \
- inc = 1 + k % (h->n_buckets - 1); last = i; \
- while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \
- if (i + inc >= h->n_buckets) i = i + inc - h->n_buckets; \
- else i += inc; \
- if (i == last) return h->n_buckets; \
- } \
- return __ac_iseither(h->flags, i)? h->n_buckets : i; \
- } else return 0; \
- } \
- static inline void kh_resize_##name(kh_##name##_t *h, khint_t new_n_buckets) \
- { \
- uint32_t *new_flags = 0; \
- khint_t j = 1; \
- { \
- khint_t t = __ac_HASH_PRIME_SIZE - 1; \
- while (__ac_prime_list[t] > new_n_buckets) --t; \
- new_n_buckets = __ac_prime_list[t+1]; \
- if (h->size >= (khint_t)(new_n_buckets * __ac_HASH_UPPER + 0.5)) j = 0; \
- else { \
- new_flags = (uint32_t*)malloc(((new_n_buckets>>4) + 1) * sizeof(uint32_t)); \
- memset(new_flags, 0xaa, ((new_n_buckets>>4) + 1) * sizeof(uint32_t)); \
- if (h->n_buckets < new_n_buckets) { \
- h->keys = (khkey_t*)realloc(h->keys, new_n_buckets * sizeof(khkey_t)); \
- if (kh_is_map) \
- h->vals = (khval_t*)realloc(h->vals, new_n_buckets * sizeof(khval_t)); \
- } \
- } \
- } \
- if (j) { \
- for (j = 0; j != h->n_buckets; ++j) { \
- if (__ac_iseither(h->flags, j) == 0) { \
- khkey_t key = h->keys[j]; \
- khval_t val; \
- if (kh_is_map) val = h->vals[j]; \
- __ac_set_isdel_true(h->flags, j); \
- while (1) { \
- khint_t inc, k, i; \
- k = __hash_func(key); \
- i = k % new_n_buckets; \
- inc = 1 + k % (new_n_buckets - 1); \
- while (!__ac_isempty(new_flags, i)) { \
- if (i + inc >= new_n_buckets) i = i + inc - new_n_buckets; \
- else i += inc; \
- } \
- __ac_set_isempty_false(new_flags, i); \
- if (i < h->n_buckets && __ac_iseither(h->flags, i) == 0) { \
- { khkey_t tmp = h->keys[i]; h->keys[i] = key; key = tmp; } \
- if (kh_is_map) { khval_t tmp = h->vals[i]; h->vals[i] = val; val = tmp; } \
- __ac_set_isdel_true(h->flags, i); \
- } else { \
- h->keys[i] = key; \
- if (kh_is_map) h->vals[i] = val; \
- break; \
- } \
- } \
- } \
- } \
- if (h->n_buckets > new_n_buckets) { \
- h->keys = (khkey_t*)realloc(h->keys, new_n_buckets * sizeof(khkey_t)); \
- if (kh_is_map) \
- h->vals = (khval_t*)realloc(h->vals, new_n_buckets * sizeof(khval_t)); \
- } \
- free(h->flags); \
- h->flags = new_flags; \
- h->n_buckets = new_n_buckets; \
- h->n_occupied = h->size; \
- h->upper_bound = (khint_t)(h->n_buckets * __ac_HASH_UPPER + 0.5); \
- } \
- } \
- static inline khint_t kh_put_##name(kh_##name##_t *h, khkey_t key, int *ret) \
- { \
- khint_t x; \
- if (h->n_occupied >= h->upper_bound) { \
- if (h->n_buckets > (h->size<<1)) kh_resize_##name(h, h->n_buckets - 1); \
- else kh_resize_##name(h, h->n_buckets + 1); \
- } \
- { \
- khint_t inc, k, i, site, last; \
- x = site = h->n_buckets; k = __hash_func(key); i = k % h->n_buckets; \
- if (__ac_isempty(h->flags, i)) x = i; \
- else { \
- inc = 1 + k % (h->n_buckets - 1); last = i; \
- while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \
- if (__ac_isdel(h->flags, i)) site = i; \
- if (i + inc >= h->n_buckets) i = i + inc - h->n_buckets; \
- else i += inc; \
- if (i == last) { x = site; break; } \
- } \
- if (x == h->n_buckets) { \
- if (__ac_isempty(h->flags, i) && site != h->n_buckets) x = site; \
- else x = i; \
- } \
- } \
- } \
- if (__ac_isempty(h->flags, x)) { \
- h->keys[x] = key; \
- __ac_set_isboth_false(h->flags, x); \
- ++h->size; ++h->n_occupied; \
- *ret = 1; \
- } else if (__ac_isdel(h->flags, x)) { \
- h->keys[x] = key; \
- __ac_set_isboth_false(h->flags, x); \
- ++h->size; \
- *ret = 2; \
- } else *ret = 0; \
- return x; \
- } \
- static inline void kh_del_##name(kh_##name##_t *h, khint_t x) \
- { \
- if (x != h->n_buckets && !__ac_iseither(h->flags, x)) { \
- __ac_set_isdel_true(h->flags, x); \
- --h->size; \
- } \
- }
-
-/* --- BEGIN OF HASH FUNCTIONS --- */
-
-/*! @function
- @abstract Integer hash function
- @param key The integer [uint32_t]
- @return The hash value [khint_t]
- */
-#define kh_int_hash_func(key) (uint32_t)(key)
-/*! @function
- @abstract Integer comparison function
- */
-#define kh_int_hash_equal(a, b) ((a) == (b))
-/*! @function
- @abstract 64-bit integer hash function
- @param key The integer [uint64_t]
- @return The hash value [khint_t]
- */
-#define kh_int64_hash_func(key) (uint32_t)((key)>>33^(key)^(key)<<11)
-/*! @function
- @abstract 64-bit integer comparison function
- */
-#define kh_int64_hash_equal(a, b) ((a) == (b))
-/*! @function
- @abstract const char* hash function
- @param s Pointer to a null terminated string
- @return The hash value
- */
-static inline khint_t __ac_X31_hash_string(const char *s)
-{
- khint_t h = *s;
- if (h) for (++s ; *s; ++s) h = (h << 5) - h + *s;
- return h;
-}
-/*! @function
- @abstract Another interface to const char* hash function
- @param key Pointer to a null terminated string [const char*]
- @return The hash value [khint_t]
- */
-#define kh_str_hash_func(key) __ac_X31_hash_string(key)
-/*! @function
- @abstract Const char* comparison function
- */
-#define kh_str_hash_equal(a, b) (strcmp(a, b) == 0)
-
-/* --- END OF HASH FUNCTIONS --- */
-
-/* Other necessary macros... */
-
-/*!
- @abstract Type of the hash table.
- @param name Name of the hash table [symbol]
- */
-#define khash_t(name) kh_##name##_t
-
-/*! @function
- @abstract Initiate a hash table.
- @param name Name of the hash table [symbol]
- @return Pointer to the hash table [khash_t(name)*]
- */
-#define kh_init(name) kh_init_##name()
-
-/*! @function
- @abstract Destroy a hash table.
- @param name Name of the hash table [symbol]
- @param h Pointer to the hash table [khash_t(name)*]
- */
-#define kh_destroy(name, h) kh_destroy_##name(h)
-
-/*! @function
- @abstract Reset a hash table without deallocating memory.
- @param name Name of the hash table [symbol]
- @param h Pointer to the hash table [khash_t(name)*]
- */
-#define kh_clear(name, h) kh_clear_##name(h)
-
-/*! @function
- @abstract Resize a hash table.
- @param name Name of the hash table [symbol]
- @param h Pointer to the hash table [khash_t(name)*]
- @param s New size [khint_t]
- */
-#define kh_resize(name, h, s) kh_resize_##name(h, s)
-
-/*! @function
- @abstract Insert a key to the hash table.
- @param name Name of the hash table [symbol]
- @param h Pointer to the hash table [khash_t(name)*]
- @param k Key [type of keys]
- @param r Extra return code: 0 if the key is present in the hash table;
- 1 if the bucket is empty (never used); 2 if the element in
- the bucket has been deleted [int*]
- @return Iterator to the inserted element [khint_t]
- */
-#define kh_put(name, h, k, r) kh_put_##name(h, k, r)
-
-/*! @function
- @abstract Retrieve a key from the hash table.
- @param name Name of the hash table [symbol]
- @param h Pointer to the hash table [khash_t(name)*]
- @param k Key [type of keys]
- @return Iterator to the found element, or kh_end(h) is the element is absent [khint_t]
- */
-#define kh_get(name, h, k) kh_get_##name(h, k)
-
-/*! @function
- @abstract Remove a key from the hash table.
- @param name Name of the hash table [symbol]
- @param h Pointer to the hash table [khash_t(name)*]
- @param k Iterator to the element to be deleted [khint_t]
- */
-#define kh_del(name, h, k) kh_del_##name(h, k)
-
-
-/*! @function
- @abstract Test whether a bucket contains data.
- @param h Pointer to the hash table [khash_t(name)*]
- @param x Iterator to the bucket [khint_t]
- @return 1 if containing data; 0 otherwise [int]
- */
-#define kh_exist(h, x) (!__ac_iseither((h)->flags, (x)))
-
-/*! @function
- @abstract Get key given an iterator
- @param h Pointer to the hash table [khash_t(name)*]
- @param x Iterator to the bucket [khint_t]
- @return Key [type of keys]
- */
-#define kh_key(h, x) ((h)->keys[x])
-
-/*! @function
- @abstract Get value given an iterator
- @param h Pointer to the hash table [khash_t(name)*]
- @param x Iterator to the bucket [khint_t]
- @return Value [type of values]
- @discussion For hash sets, calling this results in segfault.
- */
-#define kh_val(h, x) ((h)->vals[x])
-
-/*! @function
- @abstract Alias of kh_val()
- */
-#define kh_value(h, x) ((h)->vals[x])
-
-/*! @function
- @abstract Get the start iterator
- @param h Pointer to the hash table [khash_t(name)*]
- @return The start iterator [khint_t]
- */
-#define kh_begin(h) (khint_t)(0)
-
-/*! @function
- @abstract Get the end iterator
- @param h Pointer to the hash table [khash_t(name)*]
- @return The end iterator [khint_t]
- */
-#define kh_end(h) ((h)->n_buckets)
-
-/*! @function
- @abstract Get the number of elements in the hash table
- @param h Pointer to the hash table [khash_t(name)*]
- @return Number of elements in the hash table [khint_t]
- */
-#define kh_size(h) ((h)->size)
-
-/*! @function
- @abstract Get the number of buckets in the hash table
- @param h Pointer to the hash table [khash_t(name)*]
- @return Number of buckets in the hash table [khint_t]
- */
-#define kh_n_buckets(h) ((h)->n_buckets)
-
-/* More conenient interfaces */
-
-/*! @function
- @abstract Instantiate a hash set containing integer keys
- @param name Name of the hash table [symbol]
- */
-#define KHASH_SET_INIT_INT(name) \
- KHASH_INIT(name, uint32_t, char, 0, kh_int_hash_func, kh_int_hash_equal)
-
-/*! @function
- @abstract Instantiate a hash map containing integer keys
- @param name Name of the hash table [symbol]
- @param khval_t Type of values [type]
- */
-#define KHASH_MAP_INIT_INT(name, khval_t) \
- KHASH_INIT(name, uint32_t, khval_t, 1, kh_int_hash_func, kh_int_hash_equal)
-
-/*! @function
- @abstract Instantiate a hash map containing 64-bit integer keys
- @param name Name of the hash table [symbol]
- */
-#define KHASH_SET_INIT_INT64(name) \
- KHASH_INIT(name, uint64_t, char, 0, kh_int64_hash_func, kh_int64_hash_equal)
-
-/*! @function
- @abstract Instantiate a hash map containing 64-bit integer keys
- @param name Name of the hash table [symbol]
- @param khval_t Type of values [type]
- */
-#define KHASH_MAP_INIT_INT64(name, khval_t) \
- KHASH_INIT(name, uint64_t, khval_t, 1, kh_int64_hash_func, kh_int64_hash_equal)
-
-typedef const char *kh_cstr_t;
-/*! @function
- @abstract Instantiate a hash map containing const char* keys
- @param name Name of the hash table [symbol]
- */
-#define KHASH_SET_INIT_STR(name) \
- KHASH_INIT(name, kh_cstr_t, char, 0, kh_str_hash_func, kh_str_hash_equal)
-
-/*! @function
- @abstract Instantiate a hash map containing const char* keys
- @param name Name of the hash table [symbol]
- @param khval_t Type of values [type]
- */
-#define KHASH_MAP_INIT_STR(name, khval_t) \
- KHASH_INIT(name, kh_cstr_t, khval_t, 1, kh_str_hash_func, kh_str_hash_equal)
-
-#endif /* __AC_KHASH_H */
diff --git a/knetfile.c b/knetfile.c
deleted file mode 100644
index 7c96a3e..0000000
--- a/knetfile.c
+++ /dev/null
@@ -1,632 +0,0 @@
-/* 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 at 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 <ctype.h>
-#include <stdlib.h>
-#include <string.h>
-#include <errno.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>
-#endif
-
-#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;
- struct timeval tv;
- int ret;
- tv.tv_sec = 5; tv.tv_usec = 0; // 5 seconds time out
- FD_ZERO(&fds);
- FD_SET(fd, &fds);
- if (is_read) fdr = &fds;
- else fdw = &fds;
- ret = select(fd+1, fdr, fdw, 0, &tv);
-#ifndef _WIN32
- if (ret == -1) perror("select");
-#else
- if (ret == 0)
- fprintf(stderr, "select time-out\n");
- else if (ret == SOCKET_ERROR)
- fprintf(stderr, "select: %d\n", WSAGetLastError());
-#endif
- return ret;
-}
-
-#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 *int64tostr(char *buf, int64_t x)
-{
- int cnt;
- int i = 0;
- do {
- buf[i++] = '0' + x % 10;
- x /= 10;
- } while (x);
- buf[i] = 0;
- for (cnt = i, i = 0; i < cnt/2; ++i) {
- int c = buf[i]; buf[i] = buf[cnt-i-1]; buf[cnt-i-1] = c;
- }
- return buf;
-}
-
-int64_t strtoint64(const char *buf)
-{
- int64_t x;
- for (x = 0; *buf != '\0'; ++buf)
- x = x * 10 + ((int64_t) *buf - 48);
- return x;
-}
-/* In windows, the first thing is to establish the TCP connection. */
-int knet_win32_init()
-{
- 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 { \
- fprintf(stderr, "%s: %d\n", func, WSAGetLastError()); \
- 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)
-{
-#ifndef _WIN32
- unsigned char c;
-#else
- char c;
-#endif
- int n = 0;
- char *p;
- if (socket_wait(ftp->ctrl_fd, 1) <= 0) return 0;
- 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;
- ftp->response = realloc(ftp->response, ftp->max_response);
- }
- ftp->response[n++] = c;
- if (c == '\n') {
- if (n >= 4 && isdigit(ftp->response[0]) && isdigit(ftp->response[1]) && isdigit(ftp->response[2])
- && ftp->response[3] != '-') break;
- n = 0;
- continue;
- }
- }
- if (n < 2) return -1;
- ftp->response[n-2] = 0;
- return strtol(ftp->response, &p, 0);
-}
-
-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
- netwrite(ftp->ctrl_fd, cmd, strlen(cmd));
- return is_get? kftp_get_response(ftp) : 0;
-}
-
-static int kftp_pasv_prep(knetFile *ftp)
-{
- char *p;
- int v[6];
- kftp_send_cmd(ftp, "PASV\r\n", 1);
- for (p = ftp->response; *p && *p != '('; ++p);
- if (*p != '(') return -1;
- ++p;
- sscanf(p, "%d,%d,%d,%d,%d,%d", &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
- memcpy(ftp->pasv_ip, v, 4 * sizeof(int));
- ftp->pasv_port = (v[4]<<8&0xff00) + v[5];
- return 0;
-}
-
-
-static int kftp_pasv_connect(knetFile *ftp)
-{
- 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;
- }
- 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);
- ftp->fd = socket_connect(host, port);
- if (ftp->fd == -1) return -1;
- return 0;
-}
-
-int kftp_connect(knetFile *ftp)
-{
- 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)
-{
- if (ftp->ctrl_fd != -1) {
- netclose(ftp->ctrl_fd);
- ftp->ctrl_fd = -1;
- }
- netclose(ftp->fd);
- ftp->fd = -1;
- return kftp_connect(ftp);
-}
-
-// initialize ->type, ->host, ->retr and ->size
-knetFile *kftp_parse_url(const char *fn, const char *mode)
-{
- knetFile *fp;
- char *p;
- int l;
- if (strstr(fn, "ftp://") != fn) return 0;
- for (p = (char*)fn + 6; *p && *p != '/'; ++p);
- if (*p != '/') return 0;
- l = p - fn - 6;
- 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->retr = calloc(strlen(p) + 8, 1);
- sprintf(fp->retr, "RETR %s\r\n", p);
- fp->size_cmd = calloc(strlen(p) + 8, 1);
- sprintf(fp->size_cmd, "SIZE %s\r\n", p);
- fp->seek_offset = 0;
- return fp;
-}
-// place ->fd at offset off
-int kftp_connect_file(knetFile *fp)
-{
- int ret;
- long long file_size;
- if (fp->fd != -1) {
- netclose(fp->fd);
- if (fp->no_reconnect) kftp_get_response(fp);
- }
- kftp_pasv_prep(fp);
- kftp_send_cmd(fp, fp->size_cmd, 1);
-#ifndef _WIN32
- if ( sscanf(fp->response,"%*d %lld", &file_size) != 1 )
- {
- fprintf(stderr,"[kftp_connect_file] %s\n", fp->response);
- return -1;
- }
-#else
- const char *p = fp->response;
- while (*p != ' ') ++p;
- while (*p < '0' || *p > '9') ++p;
- file_size = strtoint64(p);
-#endif
- fp->file_size = file_size;
- if (fp->offset>=0) {
- char tmp[32];
-#ifndef _WIN32
- sprintf(tmp, "REST %lld\r\n", (long long)fp->offset);
-#else
- strcpy(tmp, "REST ");
- int64tostr(tmp + 5, fp->offset);
- strcat(tmp, "\r\n");
-#endif
- kftp_send_cmd(fp, tmp, 1);
- }
- kftp_send_cmd(fp, fp->retr, 0);
- kftp_pasv_connect(fp);
- ret = kftp_get_response(fp);
- if (ret != 150) {
- fprintf(stderr, "[kftp_connect_file] %s\n", fp->response);
- netclose(fp->fd);
- fp->fd = -1;
- return -1;
- }
- fp->is_ready = 1;
- 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 = 0;
- 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);
- 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>0) { // 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;
- if (mode[0] != 'r') {
- fprintf(stderr, "[kftp_open] only mode \"r\" is supported.\n");
- return 0;
- }
- if (strstr(fn, "ftp://") == fn) {
- fp = kftp_parse_url(fn, mode);
- if (fp == 0) return 0;
- if (kftp_connect(fp) == -1) {
- knet_close(fp);
- return 0;
- }
- kftp_connect_file(fp);
- } 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);
-#endif
- if (fd == -1) {
- perror("open");
- return 0;
- }
- 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;
-}
-
-knetFile *knet_dopen(int fd, const char *mode)
-{
- knetFile *fp = (knetFile*)calloc(1, sizeof(knetFile));
- fp->type = KNF_TYPE_LOCAL;
- fp->fd = fd;
- return fp;
-}
-
-off_t knet_read(knetFile *fp, void *buf, off_t len)
-{
- off_t l = 0;
- 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);
- }
- } 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) {
- curr = read(fp->fd, buf + l, rest);
- if (curr == 0) break;
- l += curr; rest -= curr;
- }
- } else l = my_netread(fp->fd, buf, len);
- fp->offset += l;
- return l;
-}
-
-off_t knet_seek(knetFile *fp, int64_t off, int whence)
-{
- if (whence == SEEK_SET && off == fp->offset) return 0;
- if (fp->type == KNF_TYPE_LOCAL) {
- /* 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) {
- // Be silent, it is OK for knet_seek to fail when the file is streamed
- // fprintf(stderr,"[knet_seek] %s\n", strerror(errno));
- return -1;
- }
- fp->offset = offset;
- return 0;
- }
- else if (fp->type == KNF_TYPE_FTP)
- {
- if (whence==SEEK_CUR)
- fp->offset += off;
- else if (whence==SEEK_SET)
- fp->offset = off;
- else if ( whence==SEEK_END)
- fp->offset = fp->file_size+off;
- fp->is_ready = 0;
- return 0;
- }
- else if (fp->type == KNF_TYPE_HTTP)
- {
- if (whence == SEEK_END) { // FIXME: can we allow SEEK_END in future?
- fprintf(stderr, "[knet_seek] SEEK_END is not supported for HTTP. Offset is unchanged.\n");
- errno = ESPIPE;
- return -1;
- }
- if (whence==SEEK_CUR)
- fp->offset += off;
- else if (whence==SEEK_SET)
- fp->offset = off;
- fp->is_ready = 0;
- return fp->offset;
- }
- errno = EINVAL;
- fprintf(stderr,"[knet_seek] %s\n", strerror(errno));
- return -1;
-}
-
-int knet_close(knetFile *fp)
-{
- if (fp == 0) return 0;
- 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); free(fp->size_cmd); // FTP specific
- free(fp->path); free(fp->http_host); // HTTP specific
- free(fp);
- return 0;
-}
-
-#ifdef KNETFILE_MAIN
-int main(void)
-{
- char *buf;
- knetFile *fp;
- 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);
- free(buf);
- return 0;
-}
-#endif
diff --git a/knetfile.h b/knetfile.h
deleted file mode 100644
index 0a0e66f..0000000
--- a/knetfile.h
+++ /dev/null
@@ -1,75 +0,0 @@
-#ifndef KNETFILE_H
-#define KNETFILE_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 <winsock2.h>
-#define netread(fd, ptr, len) recv(fd, ptr, len, 0)
-#define netwrite(fd, ptr, len) send(fd, ptr, len, 0)
-#define netclose(fd) closesocket(fd)
-#endif
-
-// FIXME: currently I/O is unbuffered
-
-#define KNF_TYPE_LOCAL 1
-#define KNF_TYPE_FTP 2
-#define KNF_TYPE_HTTP 3
-
-typedef struct knetFile_s {
- int type, fd;
- int64_t offset;
- 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, *size_cmd;
- int64_t seek_offset; // for lazy seek
- int64_t file_size;
-
- // the following are for HTTP only
- char *path, *http_host;
-} knetFile;
-
-#define knet_tell(fp) ((fp)->offset)
-#define knet_fileno(fp) ((fp)->fd)
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#ifdef _WIN32
- int knet_win32_init();
- void knet_win32_destroy();
-#endif
-
- knetFile *knet_open(const char *fn, const char *mode);
-
- /*
- This only works with local files.
- */
- knetFile *knet_dopen(int fd, const char *mode);
-
- /*
- If ->is_ready==0, this routine updates ->fd; otherwise, it simply
- reads from ->fd.
- */
- off_t knet_read(knetFile *fp, void *buf, off_t len);
-
- /*
- This routine only sets ->offset and ->is_ready=0. It does not
- communicate with the FTP server.
- */
- off_t knet_seek(knetFile *fp, int64_t off, int whence);
- int knet_close(knetFile *fp);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
diff --git a/kseq.h b/kseq.h
deleted file mode 100644
index 82face0..0000000
--- a/kseq.h
+++ /dev/null
@@ -1,227 +0,0 @@
-/* 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 at sanger.ac.uk> */
-
-/*
- 2009-07-16 (lh3): in kstream_t, change "char*" to "unsigned char*"
- */
-
-/* Last Modified: 12APR2009 */
-
-#ifndef AC_KSEQ_H
-#define AC_KSEQ_H
-
-#include <ctype.h>
-#include <string.h>
-#include <stdlib.h>
-
-#define KS_SEP_SPACE 0 // isspace(): \t, \n, \v, \f, \r
-#define KS_SEP_TAB 1 // isspace() && !' '
-#define KS_SEP_MAX 1
-
-#define __KS_TYPE(type_t) \
- typedef struct __kstream_t { \
- unsigned char *buf; \
- int begin, end, is_eof; \
- type_t f; \
- } kstream_t;
-
-#define ks_eof(ks) ((ks)->is_eof && (ks)->begin >= (ks)->end)
-#define ks_rewind(ks) ((ks)->is_eof = (ks)->begin = (ks)->end = 0)
-
-#define __KS_BASIC(type_t, __bufsize) \
- static inline kstream_t *ks_init(type_t f) \
- { \
- kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t)); \
- ks->f = f; \
- ks->buf = malloc(__bufsize); \
- return ks; \
- } \
- static inline void ks_destroy(kstream_t *ks) \
- { \
- if (ks) { \
- free(ks->buf); \
- free(ks); \
- } \
- }
-
-#define __KS_GETC(__read, __bufsize) \
- static inline int ks_getc(kstream_t *ks) \
- { \
- if (ks->is_eof && ks->begin >= ks->end) return -1; \
- if (ks->begin >= ks->end) { \
- ks->begin = 0; \
- ks->end = __read(ks->f, ks->buf, __bufsize); \
- if (ks->end < __bufsize) ks->is_eof = 1; \
- if (ks->end == 0) return -1; \
- } \
- return (int)ks->buf[ks->begin++]; \
- }
-
-#ifndef KSTRING_T
-#define KSTRING_T kstring_t
-typedef struct __kstring_t {
- size_t l, m;
- char *s;
-} kstring_t;
-#endif
-
-#ifndef kroundup32
-#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
-#endif
-
-#define __KS_GETUNTIL(__read, __bufsize) \
- static int ks_getuntil(kstream_t *ks, int delimiter, kstring_t *str, int *dret) \
- { \
- if (dret) *dret = 0; \
- str->l = 0; \
- if (ks->begin >= ks->end && ks->is_eof) return -1; \
- for (;;) { \
- int i; \
- if (ks->begin >= ks->end) { \
- if (!ks->is_eof) { \
- ks->begin = 0; \
- ks->end = __read(ks->f, ks->buf, __bufsize); \
- if (ks->end < __bufsize) ks->is_eof = 1; \
- if (ks->end == 0) break; \
- } else break; \
- } \
- if (delimiter > KS_SEP_MAX) { \
- for (i = ks->begin; i < ks->end; ++i) \
- if (ks->buf[i] == delimiter) break; \
- } else if (delimiter == KS_SEP_SPACE) { \
- for (i = ks->begin; i < ks->end; ++i) \
- if (isspace(ks->buf[i])) break; \
- } else if (delimiter == KS_SEP_TAB) { \
- for (i = ks->begin; i < ks->end; ++i) \
- if (isspace(ks->buf[i]) && ks->buf[i] != ' ') break; \
- } else i = 0; /* never come to here! */ \
- if (str->m - str->l < i - ks->begin + 1) { \
- str->m = str->l + (i - ks->begin) + 1; \
- kroundup32(str->m); \
- str->s = (char*)realloc(str->s, str->m); \
- } \
- memcpy(str->s + str->l, ks->buf + ks->begin, i - ks->begin); \
- str->l = str->l + (i - ks->begin); \
- ks->begin = i + 1; \
- if (i < ks->end) { \
- if (dret) *dret = ks->buf[i]; \
- break; \
- } \
- } \
- if (str->l == 0) { \
- str->m = 1; \
- str->s = (char*)calloc(1, 1); \
- } \
- str->s[str->l] = '\0'; \
- return str->l; \
- }
-
-#define KSTREAM_INIT(type_t, __read, __bufsize) \
- __KS_TYPE(type_t) \
- __KS_BASIC(type_t, __bufsize) \
- __KS_GETC(__read, __bufsize) \
- __KS_GETUNTIL(__read, __bufsize)
-
-#define __KSEQ_BASIC(type_t) \
- static inline kseq_t *kseq_init(type_t fd) \
- { \
- kseq_t *s = (kseq_t*)calloc(1, sizeof(kseq_t)); \
- s->f = ks_init(fd); \
- return s; \
- } \
- static inline void kseq_rewind(kseq_t *ks) \
- { \
- ks->last_char = 0; \
- ks->f->is_eof = ks->f->begin = ks->f->end = 0; \
- } \
- static inline void kseq_destroy(kseq_t *ks) \
- { \
- if (!ks) return; \
- free(ks->name.s); free(ks->comment.s); free(ks->seq.s); free(ks->qual.s); \
- ks_destroy(ks->f); \
- free(ks); \
- }
-
-/* Return value:
- >=0 length of the sequence (normal)
- -1 end-of-file
- -2 truncated quality string
- */
-#define __KSEQ_READ \
- static int kseq_read(kseq_t *seq) \
- { \
- int c; \
- kstream_t *ks = seq->f; \
- if (seq->last_char == 0) { /* then jump to the next header line */ \
- while ((c = ks_getc(ks)) != -1 && c != '>' && c != '@'); \
- if (c == -1) return -1; /* end of file */ \
- seq->last_char = c; \
- } /* the first header char has been read */ \
- seq->comment.l = seq->seq.l = seq->qual.l = 0; \
- if (ks_getuntil(ks, 0, &seq->name, &c) < 0) return -1; \
- if (c != '\n') ks_getuntil(ks, '\n', &seq->comment, 0); \
- while ((c = ks_getc(ks)) != -1 && c != '>' && c != '+' && c != '@') { \
- if (isgraph(c)) { /* printable non-space character */ \
- if (seq->seq.l + 1 >= seq->seq.m) { /* double the memory */ \
- seq->seq.m = seq->seq.l + 2; \
- kroundup32(seq->seq.m); /* rounded to next closest 2^k */ \
- seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \
- } \
- seq->seq.s[seq->seq.l++] = (char)c; \
- } \
- } \
- if (c == '>' || c == '@') seq->last_char = c; /* the first header char has been read */ \
- seq->seq.s[seq->seq.l] = 0; /* null terminated string */ \
- if (c != '+') return seq->seq.l; /* FASTA */ \
- if (seq->qual.m < seq->seq.m) { /* allocate enough memory */ \
- seq->qual.m = seq->seq.m; \
- seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m); \
- } \
- while ((c = ks_getc(ks)) != -1 && c != '\n'); /* skip the rest of '+' line */ \
- if (c == -1) return -2; /* we should not stop here */ \
- while ((c = ks_getc(ks)) != -1 && seq->qual.l < seq->seq.l) \
- if (c >= 33 && c <= 127) seq->qual.s[seq->qual.l++] = (unsigned char)c; \
- seq->qual.s[seq->qual.l] = 0; /* null terminated string */ \
- seq->last_char = 0; /* we have not come to the next header line */ \
- if (seq->seq.l != seq->qual.l) return -2; /* qual string is shorter than seq string */ \
- return seq->seq.l; \
- }
-
-#define __KSEQ_TYPE(type_t) \
- typedef struct { \
- kstring_t name, comment, seq, qual; \
- int last_char; \
- kstream_t *f; \
- } kseq_t;
-
-#define KSEQ_INIT(type_t, __read) \
- KSTREAM_INIT(type_t, __read, 4096) \
- __KSEQ_TYPE(type_t) \
- __KSEQ_BASIC(type_t) \
- __KSEQ_READ
-
-#endif
diff --git a/ksort.h b/ksort.h
deleted file mode 100644
index 16a03fd..0000000
--- a/ksort.h
+++ /dev/null
@@ -1,271 +0,0 @@
-/* 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 at sanger.ac.uk> */
-
-/*
- 2008-11-16 (0.1.4):
-
- * Fixed a bug in introsort() that happens in rare cases.
-
- 2008-11-05 (0.1.3):
-
- * Fixed a bug in introsort() for complex comparisons.
-
- * Fixed a bug in mergesort(). The previous version is not stable.
-
- 2008-09-15 (0.1.2):
-
- * Accelerated introsort. On my Mac (not on another Linux machine),
- my implementation is as fast as std::sort on random input.
-
- * Added combsort and in introsort, switch to combsort if the
- recursion is too deep.
-
- 2008-09-13 (0.1.1):
-
- * Added k-small algorithm
-
- 2008-09-05 (0.1.0):
-
- * Initial version
-
-*/
-
-#ifndef AC_KSORT_H
-#define AC_KSORT_H
-
-#include <stdlib.h>
-#include <string.h>
-
-typedef struct {
- void *left, *right;
- int depth;
-} ks_isort_stack_t;
-
-#define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
-
-#define KSORT_INIT(name, type_t, __sort_lt) \
- void ks_mergesort_##name(size_t n, type_t array[], type_t temp[]) \
- { \
- type_t *a2[2], *a, *b; \
- int curr, shift; \
- \
- a2[0] = array; \
- a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n); \
- for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) { \
- a = a2[curr]; b = a2[1-curr]; \
- if (shift == 0) { \
- type_t *p = b, *i, *eb = a + n; \
- for (i = a; i < eb; i += 2) { \
- if (i == eb - 1) *p++ = *i; \
- else { \
- if (__sort_lt(*(i+1), *i)) { \
- *p++ = *(i+1); *p++ = *i; \
- } else { \
- *p++ = *i; *p++ = *(i+1); \
- } \
- } \
- } \
- } else { \
- size_t i, step = 1ul<<shift; \
- for (i = 0; i < n; i += step<<1) { \
- type_t *p, *j, *k, *ea, *eb; \
- if (n < i + step) { \
- ea = a + n; eb = a; \
- } else { \
- ea = a + i + step; \
- eb = a + (n < i + (step<<1)? n : i + (step<<1)); \
- } \
- j = a + i; k = a + i + step; p = b + i; \
- while (j < ea && k < eb) { \
- if (__sort_lt(*k, *j)) *p++ = *k++; \
- else *p++ = *j++; \
- } \
- while (j < ea) *p++ = *j++; \
- while (k < eb) *p++ = *k++; \
- } \
- } \
- curr = 1 - curr; \
- } \
- if (curr == 1) { \
- type_t *p = a2[0], *i = a2[1], *eb = array + n; \
- for (; p < eb; ++i) *p++ = *i; \
- } \
- if (temp == 0) free(a2[1]); \
- } \
- void ks_heapadjust_##name(size_t i, size_t n, type_t l[]) \
- { \
- size_t k = i; \
- type_t tmp = l[i]; \
- while ((k = (k << 1) + 1) < n) { \
- if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k; \
- if (__sort_lt(l[k], tmp)) break; \
- l[i] = l[k]; i = k; \
- } \
- l[i] = tmp; \
- } \
- void ks_heapmake_##name(size_t lsize, type_t l[]) \
- { \
- size_t i; \
- for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \
- ks_heapadjust_##name(i, lsize, l); \
- } \
- void ks_heapsort_##name(size_t lsize, type_t l[]) \
- { \
- size_t i; \
- for (i = lsize - 1; i > 0; --i) { \
- type_t tmp; \
- tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \
- } \
- } \
- inline void __ks_insertsort_##name(type_t *s, type_t *t) \
- { \
- type_t *i, *j, swap_tmp; \
- for (i = s + 1; i < t; ++i) \
- for (j = i; j > s && __sort_lt(*j, *(j-1)); --j) { \
- swap_tmp = *j; *j = *(j-1); *(j-1) = swap_tmp; \
- } \
- } \
- void ks_combsort_##name(size_t n, type_t a[]) \
- { \
- const double shrink_factor = 1.2473309501039786540366528676643; \
- int do_swap; \
- size_t gap = n; \
- type_t tmp, *i, *j; \
- do { \
- if (gap > 2) { \
- gap = (size_t)(gap / shrink_factor); \
- if (gap == 9 || gap == 10) gap = 11; \
- } \
- do_swap = 0; \
- for (i = a; i < a + n - gap; ++i) { \
- j = i + gap; \
- if (__sort_lt(*j, *i)) { \
- tmp = *i; *i = *j; *j = tmp; \
- do_swap = 1; \
- } \
- } \
- } while (do_swap || gap > 2); \
- if (gap != 1) __ks_insertsort_##name(a, a + n); \
- } \
- void ks_introsort_##name(size_t n, type_t a[]) \
- { \
- int d; \
- ks_isort_stack_t *top, *stack; \
- type_t rp, swap_tmp; \
- type_t *s, *t, *i, *j, *k; \
- \
- if (n < 1) return; \
- else if (n == 2) { \
- if (__sort_lt(a[1], a[0])) { swap_tmp = a[0]; a[0] = a[1]; a[1] = swap_tmp; } \
- return; \
- } \
- for (d = 2; 1ul<<d < n; ++d); \
- stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \
- top = stack; s = a; t = a + (n-1); d <<= 1; \
- while (1) { \
- if (s < t) { \
- if (--d == 0) { \
- ks_combsort_##name(t - s + 1, s); \
- t = s; \
- continue; \
- } \
- i = s; j = t; k = i + ((j-i)>>1) + 1; \
- if (__sort_lt(*k, *i)) { \
- if (__sort_lt(*k, *j)) k = j; \
- } else k = __sort_lt(*j, *i)? i : j; \
- rp = *k; \
- if (k != t) { swap_tmp = *k; *k = *t; *t = swap_tmp; } \
- for (;;) { \
- do ++i; while (__sort_lt(*i, rp)); \
- do --j; while (i <= j && __sort_lt(rp, *j)); \
- if (j <= i) break; \
- swap_tmp = *i; *i = *j; *j = swap_tmp; \
- } \
- swap_tmp = *i; *i = *t; *t = swap_tmp; \
- if (i-s > t-i) { \
- if (i-s > 16) { top->left = s; top->right = i-1; top->depth = d; ++top; } \
- s = t-i > 16? i+1 : t; \
- } else { \
- if (t-i > 16) { top->left = i+1; top->right = t; top->depth = d; ++top; } \
- t = i-s > 16? i-1 : s; \
- } \
- } else { \
- if (top == stack) { \
- free(stack); \
- __ks_insertsort_##name(a, a+n); \
- return; \
- } else { --top; s = (type_t*)top->left; t = (type_t*)top->right; d = top->depth; } \
- } \
- } \
- } \
- /* This function is adapted from: http://ndevilla.free.fr/median/ */ \
- /* 0 <= kk < n */ \
- type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk) \
- { \
- type_t *low, *high, *k, *ll, *hh, *mid; \
- low = arr; high = arr + n - 1; k = arr + kk; \
- for (;;) { \
- if (high <= low) return *k; \
- if (high == low + 1) { \
- if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
- return *k; \
- } \
- mid = low + (high - low) / 2; \
- if (__sort_lt(*high, *mid)) KSORT_SWAP(type_t, *mid, *high); \
- if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
- if (__sort_lt(*low, *mid)) KSORT_SWAP(type_t, *mid, *low); \
- KSORT_SWAP(type_t, *mid, *(low+1)); \
- ll = low + 1; hh = high; \
- for (;;) { \
- do ++ll; while (__sort_lt(*ll, *low)); \
- do --hh; while (__sort_lt(*low, *hh)); \
- if (hh < ll) break; \
- KSORT_SWAP(type_t, *ll, *hh); \
- } \
- KSORT_SWAP(type_t, *low, *hh); \
- if (hh <= k) low = ll; \
- if (hh >= k) high = hh - 1; \
- } \
- }
-
-#define ks_mergesort(name, n, a, t) ks_mergesort_##name(n, a, t)
-#define ks_introsort(name, n, a) ks_introsort_##name(n, a)
-#define ks_combsort(name, n, a) ks_combsort_##name(n, a)
-#define ks_heapsort(name, n, a) ks_heapsort_##name(n, a)
-#define ks_heapmake(name, n, a) ks_heapmake_##name(n, a)
-#define ks_heapadjust(name, i, n, a) ks_heapadjust_##name(i, n, a)
-#define ks_ksmall(name, n, a, k) ks_ksmall_##name(n, a, k)
-
-#define ks_lt_generic(a, b) ((a) < (b))
-#define ks_lt_str(a, b) (strcmp((a), (b)) < 0)
-
-typedef const char *ksstr_t;
-
-#define KSORT_INIT_GENERIC(type_t) KSORT_INIT(type_t, type_t, ks_lt_generic)
-#define KSORT_INIT_STR KSORT_INIT(str, ksstr_t, ks_lt_str)
-
-#endif
diff --git a/kstring.c b/kstring.c
deleted file mode 100644
index e0203fa..0000000
--- a/kstring.c
+++ /dev/null
@@ -1,165 +0,0 @@
-#include <stdarg.h>
-#include <stdio.h>
-#include <ctype.h>
-#include <string.h>
-#include <stdint.h>
-#include "kstring.h"
-
-int ksprintf(kstring_t *s, const char *fmt, ...)
-{
- va_list ap;
- int l;
- va_start(ap, fmt);
- l = vsnprintf(s->s + s->l, s->m - s->l, fmt, ap); // This line does not work with glibc 2.0. See `man snprintf'.
- va_end(ap);
- if (l + 1 > s->m - s->l) {
- s->m = s->l + l + 2;
- kroundup32(s->m);
- s->s = (char*)realloc(s->s, s->m);
- va_start(ap, fmt);
- l = vsnprintf(s->s + s->l, s->m - s->l, fmt, ap);
- }
- va_end(ap);
- s->l += l;
- return l;
-}
-
-// s MUST BE a null terminated string; l = strlen(s)
-int ksplit_core(char *s, int delimiter, int *_max, int **_offsets)
-{
- int i, n, max, last_char, last_start, *offsets, l;
- n = 0; max = *_max; offsets = *_offsets;
- l = strlen(s);
-
-#define __ksplit_aux do { \
- if (_offsets) { \
- s[i] = 0; \
- if (n == max) { \
- max = max? max<<1 : 2; \
- offsets = (int*)realloc(offsets, sizeof(int) * max); \
- } \
- offsets[n++] = last_start; \
- } else ++n; \
- } while (0)
-
- for (i = 0, last_char = last_start = 0; i <= l; ++i) {
- if (delimiter == 0) {
- if (isspace(s[i]) || s[i] == 0) {
- if (isgraph(last_char)) __ksplit_aux; // the end of a field
- } else {
- if (isspace(last_char) || last_char == 0) last_start = i;
- }
- } else {
- if (s[i] == delimiter || s[i] == 0) {
- if (last_char != 0 && last_char != delimiter) __ksplit_aux; // the end of a field
- } else {
- if (last_char == delimiter || last_char == 0) last_start = i;
- }
- }
- last_char = s[i];
- }
- *_max = max; *_offsets = offsets;
- 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()
-{
- kstring_t *s;
- int *fields, n, i;
- s = (kstring_t*)calloc(1, sizeof(kstring_t));
- // test ksprintf()
- ksprintf(s, " abcdefg: %d ", 100);
- printf("'%s'\n", s->s);
- // test ksplit()
- fields = ksplit(s, 0, &n);
- 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
diff --git a/kstring.h b/kstring.h
deleted file mode 100644
index f4e5a99..0000000
--- a/kstring.h
+++ /dev/null
@@ -1,68 +0,0 @@
-#ifndef KSTRING_H
-#define KSTRING_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))
-#endif
-
-#ifndef KSTRING_T
-#define KSTRING_T kstring_t
-typedef struct __kstring_t {
- size_t l, m;
- char *s;
-} kstring_t;
-#endif
-
-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) {
- s->m = s->l + l + 2;
- kroundup32(s->m);
- s->s = (char*)realloc(s->s, s->m);
- }
- strncpy(s->s + s->l, p, l);
- s->l += l;
- s->s[s->l] = 0;
- return l;
-}
-
-static inline int kputs(const char *p, kstring_t *s)
-{
- return kputsn(p, strlen(p), s);
-}
-
-static inline int kputc(int c, kstring_t *s)
-{
- if (s->l + 1 >= s->m) {
- s->m = s->l + 2;
- kroundup32(s->m);
- s->s = (char*)realloc(s->s, s->m);
- }
- s->s[s->l++] = c;
- s->s[s->l] = 0;
- return c;
-}
-
-static inline int *ksplit(kstring_t *s, int delimiter, int *n)
-{
- int max = 0, *offsets = 0;
- *n = ksplit_core(s->s, delimiter, &max, &offsets);
- return offsets;
-}
-
-#endif
diff --git a/main.c b/main.c
deleted file mode 100644
index 364abe5..0000000
--- a/main.c
+++ /dev/null
@@ -1,290 +0,0 @@
-#include <string.h>
-#include <unistd.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <sys/stat.h>
-#include <errno.h>
-#include "bgzf.h"
-#include "tabix.h"
-
-#define PACKAGE_VERSION "0.2.5 (r964)"
-
-#define error(...) { fprintf(stderr,__VA_ARGS__); return -1; }
-
-int reheader_file(const char *header, const char *file, int meta)
-{
- BGZF *fp = bgzf_open(file,"r");
- if (bgzf_read_block(fp) != 0 || !fp->block_length)
- return -1;
-
- char *buffer = fp->uncompressed_block;
- int skip_until = 0;
-
- if ( buffer[0]==meta )
- {
- skip_until = 1;
-
- // Skip the header
- while (1)
- {
- if ( buffer[skip_until]=='\n' )
- {
- skip_until++;
- if ( skip_until>=fp->block_length )
- {
- if (bgzf_read_block(fp) != 0 || !fp->block_length)
- error("no body?\n");
- skip_until = 0;
- }
- // The header has finished
- if ( buffer[skip_until]!=meta ) break;
- }
- skip_until++;
- if ( skip_until>=fp->block_length )
- {
- if (bgzf_read_block(fp) != 0 || !fp->block_length)
- error("no body?\n");
- skip_until = 0;
- }
- }
- }
-
- FILE *fh = fopen(header,"r");
- if ( !fh )
- error("%s: %s", header,strerror(errno));
- int page_size = getpagesize();
- char *buf = valloc(page_size);
- BGZF *bgzf_out = bgzf_fdopen(fileno(stdout), "w");
- ssize_t nread;
- while ( (nread=fread(buf,1,page_size-1,fh))>0 )
- {
- if ( nread<page_size-1 && buf[nread-1]!='\n' )
- buf[nread++] = '\n';
- if (bgzf_write(bgzf_out, buf, nread) < 0) error("Error: %s\n",bgzf_out->error);
- }
- fclose(fh);
-
- if ( fp->block_length - skip_until > 0 )
- {
- if (bgzf_write(bgzf_out, buffer+skip_until, fp->block_length-skip_until) < 0)
- error("Error: %s\n",fp->error);
- }
- if (bgzf_flush(bgzf_out) < 0)
- error("Error: %s\n",bgzf_out->error);
-
- while (1)
- {
-#ifdef _USE_KNETFILE
- nread = knet_read(fp->x.fpr, buf, page_size);
-#else
- nread = fread(buf, 1, page_size, fp->file);
-#endif
- if ( nread<=0 )
- break;
-
-#ifdef _USE_KNETFILE
- int count = fwrite(buf, 1, nread, bgzf_out->x.fpw);
-#else
- int count = fwrite(buf, 1, nread, bgzf_out->file);
-#endif
- if (count != nread)
- error("Write failed, wrote %d instead of %d bytes.\n", count,(int)nread);
- }
-
- if (bgzf_close(bgzf_out) < 0)
- error("Error: %s\n",bgzf_out->error);
-
- return 0;
-}
-
-
-int main(int argc, char *argv[])
-{
- int c, skip = -1, meta = -1, list_chrms = 0, force = 0, print_header = 0, bed_reg = 0;
- ti_conf_t conf = ti_conf_gff;
- const char *reheader = NULL;
- while ((c = getopt(argc, argv, "p:s:b:e:0S:c:lhfBr:")) >= 0) {
- switch (c) {
- case 'B': bed_reg = 1; break;
- case '0': conf.preset |= TI_FLAG_UCSC; break;
- case 'S': skip = atoi(optarg); break;
- case 'c': meta = optarg[0]; break;
- case 'p':
- if (strcmp(optarg, "gff") == 0) conf = ti_conf_gff;
- else if (strcmp(optarg, "bed") == 0) conf = ti_conf_bed;
- else if (strcmp(optarg, "sam") == 0) conf = ti_conf_sam;
- else if (strcmp(optarg, "vcf") == 0 || strcmp(optarg, "vcf4") == 0) conf = ti_conf_vcf;
- else if (strcmp(optarg, "psltbl") == 0) conf = ti_conf_psltbl;
- else {
- fprintf(stderr, "[main] unrecognized preset '%s'\n", optarg);
- return 1;
- }
- break;
- case 's': conf.sc = atoi(optarg); break;
- case 'b': conf.bc = atoi(optarg); break;
- case 'e': conf.ec = atoi(optarg); break;
- case 'l': list_chrms = 1; break;
- case 'h': print_header = 1; break;
- case 'f': force = 1; break;
- case 'r': reheader = optarg; break;
- }
- }
- if (skip >= 0) conf.line_skip = skip;
- if (meta >= 0) conf.meta_char = meta;
- if (optind == argc) {
- fprintf(stderr, "\n");
- fprintf(stderr, "Program: tabix (TAB-delimited file InderXer)\n");
- fprintf(stderr, "Version: %s\n\n", PACKAGE_VERSION);
- fprintf(stderr, "Usage: tabix <in.tab.bgz> [region1 [region2 [...]]]\n\n");
- fprintf(stderr, "Options: -p STR preset: gff, bed, sam, vcf, psltbl [gff]\n");
- fprintf(stderr, " -s INT sequence name column [1]\n");
- fprintf(stderr, " -b INT start column [4]\n");
- fprintf(stderr, " -e INT end column; can be identical to '-b' [5]\n");
- fprintf(stderr, " -S INT skip first INT lines [0]\n");
- fprintf(stderr, " -c CHAR symbol for comment/meta lines [#]\n");
- fprintf(stderr, " -r FILE replace the header with the content of FILE [null]\n");
- fprintf(stderr, " -B region1 is a BED file (entire file will be read)\n");
- fprintf(stderr, " -0 zero-based coordinate\n");
- fprintf(stderr, " -h print the header lines\n");
- fprintf(stderr, " -l list chromosome names\n");
- fprintf(stderr, " -f force to overwrite the index\n");
- fprintf(stderr, "\n");
- return 1;
- }
- if (list_chrms) {
- ti_index_t *idx;
- int i, n;
- const char **names;
- idx = ti_index_load(argv[optind]);
- if (idx == 0) {
- fprintf(stderr, "[main] fail to load the index file.\n");
- return 1;
- }
- names = ti_seqname(idx, &n);
- for (i = 0; i < n; ++i) printf("%s\n", names[i]);
- free(names);
- ti_index_destroy(idx);
- return 0;
- }
- if (reheader)
- return reheader_file(reheader,argv[optind],conf.meta_char);
-
- struct stat stat_tbi,stat_vcf;
- char *fnidx = calloc(strlen(argv[optind]) + 5, 1);
- strcat(strcpy(fnidx, argv[optind]), ".tbi");
-
- if (optind + 1 == argc) {
- if (force == 0) {
- if (stat(fnidx, &stat_tbi) == 0)
- {
- // Before complaining, check if the VCF file isn't newer. This is a common source of errors,
- // people tend not to notice that tabix failed
- stat(argv[optind], &stat_vcf);
- if ( stat_vcf.st_mtime <= stat_tbi.st_mtime )
- {
- fprintf(stderr, "[tabix] the index file exists. Please use '-f' to overwrite.\n");
- free(fnidx);
- return 1;
- }
- }
- }
- if ( bgzf_check_bgzf(argv[optind])!=1 )
- {
- fprintf(stderr,"[tabix] was bgzip used to compress this file? %s\n", argv[optind]);
- free(fnidx);
- return 1;
- }
- return ti_index_build(argv[optind], &conf);
- }
- { // retrieve
- tabix_t *t;
- // Common source of errors: new VCF is used with an old index
- stat(fnidx, &stat_tbi);
- stat(argv[optind], &stat_vcf);
- if ( force==0 && stat_vcf.st_mtime > stat_tbi.st_mtime )
- {
- fprintf(stderr, "[tabix] the index file is older than the vcf file. Please use '-f' to overwrite or reindex.\n");
- free(fnidx);
- return 1;
- }
- free(fnidx);
-
- if ((t = ti_open(argv[optind], 0)) == 0) {
- fprintf(stderr, "[main] fail to open the data file.\n");
- return 1;
- }
- if (strcmp(argv[optind+1], ".") == 0) { // retrieve all
- ti_iter_t iter;
- const char *s;
- int len;
- iter = ti_query(t, 0, 0, 0);
- while ((s = ti_read(t, iter, &len)) != 0) {
- fputs(s, stdout); fputc('\n', stdout);
- }
- ti_iter_destroy(iter);
- } else { // retrieve from specified regions
- int i, len;
- ti_iter_t iter;
- const char *s;
- const ti_conf_t *idxconf;
-
- if (ti_lazy_index_load(t) < 0 && bed_reg == 0) {
- fprintf(stderr,"[tabix] failed to load the index file.\n");
- return 1;
- }
- idxconf = ti_get_conf(t->idx);
-
- if ( print_header )
- {
- // If requested, print the header lines here
- iter = ti_query(t, 0, 0, 0);
- while ((s = ti_read(t, iter, &len)) != 0) {
- if ((int)(*s) != idxconf->meta_char) break;
- fputs(s, stdout); fputc('\n', stdout);
- }
- ti_iter_destroy(iter);
- }
- if (bed_reg) {
- extern int bed_overlap(const void *_h, const char *chr, int beg, int end);
- extern void *bed_read(const char *fn);
- extern void bed_destroy(void *_h);
-
- const ti_conf_t *conf_ = idxconf? idxconf : &conf; // use the index file if available
- void *bed = bed_read(argv[optind+1]); // load the BED file
- ti_interval_t intv;
-
- if (bed == 0) {
- fprintf(stderr, "[main] fail to read the BED file.\n");
- return 1;
- }
- iter = ti_query(t, 0, 0, 0);
- while ((s = ti_read(t, iter, &len)) != 0) {
- int c;
- ti_get_intv(conf_, len, (char*)s, &intv);
- c = *intv.se; *intv.se = '\0';
- if (bed_overlap(bed, intv.ss, intv.beg, intv.end)) {
- *intv.se = c;
- puts(s);
- }
- *intv.se = c;
- }
- ti_iter_destroy(iter);
- bed_destroy(bed);
- } else {
- for (i = optind + 1; i < argc; ++i) {
- int tid, beg, end;
- if (ti_parse_region(t->idx, argv[i], &tid, &beg, &end) == 0) {
- iter = ti_queryi(t, tid, beg, end);
- while ((s = ti_read(t, iter, &len)) != 0) {
- fputs(s, stdout); fputc('\n', stdout);
- }
- ti_iter_destroy(iter);
- }
- // else fprintf(stderr, "[main] invalid region: unknown target name or minus interval.\n");
- }
- }
- }
- ti_close(t);
- }
- return 0;
-}
diff --git a/perl/MANIFEST b/perl/MANIFEST
deleted file mode 100644
index 877da96..0000000
--- a/perl/MANIFEST
+++ /dev/null
@@ -1,8 +0,0 @@
-MANIFEST
-typemap
-Tabix.xs
-Tabix.pm
-TabixIterator.pm
-Makefile.PL
-t/01local.t
-t/02remote.t
\ No newline at end of file
diff --git a/perl/Makefile.PL b/perl/Makefile.PL
deleted file mode 100644
index 3ea6e8d..0000000
--- a/perl/Makefile.PL
+++ /dev/null
@@ -1,8 +0,0 @@
-use ExtUtils::MakeMaker;
-WriteMakefile(
- NAME => 'Tabix',
- VERSION_FROM => 'Tabix.pm',
- LIBS => ['-lz -L.. -ltabix'],
- DEFINE => '-D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE',
- INC => '-I..',
- );
diff --git a/perl/Tabix.pm b/perl/Tabix.pm
deleted file mode 100644
index fd7165d..0000000
--- a/perl/Tabix.pm
+++ /dev/null
@@ -1,76 +0,0 @@
-package Tabix;
-
-use strict;
-use warnings;
-use Carp qw/croak/;
-
-use TabixIterator;
-
-require Exporter;
-
-our @ISA = qw/Exporter/;
-our @EXPORT = qw/tabix_open tabix_close tabix_read tabix_query tabix_getnames tabix_iter_free/;
-
-our $VERSION = '0.2.0';
-
-require XSLoader;
-XSLoader::load('Tabix', $VERSION);
-
-sub new {
- my $invocant = shift;
- my %args = @_;
- $args{-data} || croak("-data argument required");
- my $class = ref($invocant) || $invocant;
- my $self = {};
- bless($self, $class);
- $self->open($args{-data}, $args{-index});
- return $self;
-}
-
-sub open {
- my ($self, $fn, $fnidx) = @_;
- $self->close;
- $self->{_fn} = $fn;
- $self->{_fnidx} = $fnidx;
- $self->{_} = $fnidx? tabix_open($fn, $fnidx) : tabix_open($fn);
-}
-
-sub close {
- my $self = shift;
- if ($self->{_}) {
- tabix_close($self->{_});
- delete($self->{_}); delete($self->{_fn}); delete($self->{_fnidx});
- }
-}
-
-sub DESTROY {
- my $self = shift;
- $self->close;
-}
-
-sub query {
- my $self = shift;
- my $iter;
- if (@_) {
- $iter = tabix_query($self->{_}, @_);
- } else {
- $iter = tabix_query($self->{_});
- }
- my $i = TabixIterator->new;
- $i->set($iter);
- return $i;
-}
-
-sub read {
- my $self = shift;
- my $iter = shift;
- return tabix_read($self->{_}, $iter->get);
-}
-
-sub getnames {
- my $self = shift;
- return tabix_getnames($self->{_});
-}
-
-1;
-__END__
diff --git a/perl/Tabix.xs b/perl/Tabix.xs
deleted file mode 100644
index 50dabb1..0000000
--- a/perl/Tabix.xs
+++ /dev/null
@@ -1,71 +0,0 @@
-#include "EXTERN.h"
-#include "perl.h"
-#include "XSUB.h"
-
-#include <stdlib.h>
-#include "tabix.h"
-
-MODULE = Tabix PACKAGE = Tabix
-
-tabix_t*
-tabix_open(fn, fnidx=0)
- char *fn
- char *fnidx
- CODE:
- RETVAL = ti_open(fn, fnidx);
- OUTPUT:
- RETVAL
-
-void
-tabix_close(t)
- tabix_t *t
- CODE:
- ti_close(t);
-
-ti_iter_t
-tabix_query(t, seq=0, beg=0, end=0x7fffffff)
- tabix_t *t
- const char *seq
- int beg
- int end
- PREINIT:
- CODE:
- RETVAL = ti_query(t, seq, beg, end);
- OUTPUT:
- RETVAL
-
-SV*
-tabix_read(t, iter)
- tabix_t *t
- ti_iter_t iter
- PREINIT:
- const char *s;
- int len;
- CODE:
- s = ti_read(t, iter, &len);
- if (s == 0)
- return XSRETURN_EMPTY;
- RETVAL = newSVpv(s, len);
- OUTPUT:
- RETVAL
-
-void
-tabix_getnames(t)
- tabix_t *t
- PREINIT:
- const char **names;
- int i, n;
- PPCODE:
- ti_lazy_index_load(t);
- names = ti_seqname(t->idx, &n);
- for (i = 0; i < n; ++i)
- XPUSHs(sv_2mortal(newSVpv(names[i], 0)));
- free(names);
-
-MODULE = Tabix PACKAGE = TabixIterator
-
-void
-tabix_iter_free(iter)
- ti_iter_t iter
- CODE:
- ti_iter_destroy(iter);
diff --git a/perl/TabixIterator.pm b/perl/TabixIterator.pm
deleted file mode 100644
index 335194a..0000000
--- a/perl/TabixIterator.pm
+++ /dev/null
@@ -1,41 +0,0 @@
-package TabixIterator;
-
-use strict;
-use warnings;
-use Carp qw/croak/;
-
-require Exporter;
-
-our @ISA = qw/Exporter/;
-our @EXPORT = qw/tabix_iter_free/;
-
-our $VERSION = '0.2.0';
-
-require XSLoader;
-XSLoader::load('Tabix', $VERSION);
-
-sub new {
- my $invocant = shift;
- my $class = ref($invocant) || $invocant;
- my $self = {};
- bless($self, $class);
- return $self;
-}
-
-sub set {
- my ($self, $iter) = @_;
- $self->{_} = $iter;
-}
-
-sub get {
- my $self = shift;
- return $self->{_};
-}
-
-sub DESTROY {
- my $self = shift;
- tabix_iter_free($self->{_}) if ($self->{_});
-}
-
-1;
-__END__
diff --git a/perl/t/01local.t b/perl/t/01local.t
deleted file mode 100644
index 4eb6534..0000000
--- a/perl/t/01local.t
+++ /dev/null
@@ -1,28 +0,0 @@
-#-*-Perl-*-
-use Test::More tests => 9;
-BEGIN { use_ok('Tabix') };
-
-{ # C-like low-level interface
- my $t = tabix_open("../example.gtf.gz");
- ok($t);
- my $iter = tabix_query($t, "chr1", 0, 2000);
- ok($iter);
- $_ = 0;
- ++$_ while (tabix_read($t, $iter));
- is($_, 6);
- tabix_iter_free($iter);
- @_ = tabix_getnames($t);
- is(scalar(@_), 2);
-}
-
-{ # OOP high-level interface
- my $t = Tabix->new(-data=>"../example.gtf.gz");
- ok($t);
- my $iter = $t->query("chr1", 3000, 5000);
- ok($iter);
- $_ = 0;
- ++$_ while ($t->read($iter));
- is($_, 27);
- @_ = $t->getnames;
- is($_[1], "chr2");
-}
diff --git a/perl/t/02remote.t b/perl/t/02remote.t
deleted file mode 100644
index 0668e8f..0000000
--- a/perl/t/02remote.t
+++ /dev/null
@@ -1,28 +0,0 @@
-#-*-Perl-*-
-use Test::More tests => 9;
-BEGIN { use_ok('Tabix') };
-
-{ # FTP access
- my $t = Tabix->new(-data=>"ftp://ftp.ncbi.nih.gov/1000genomes/ftp/pilot_data/release/2010_03/pilot1/CEU.SRP000031.2010_03.genotypes.vcf.gz");
- ok($t);
- my $iter = $t->query("1", 1000000, 1100000);
- ok($iter);
- $_ = 0;
- ++$_ while ($t->read($iter));
- is($_, 306);
- @_ = $t->getnames;
- is(scalar(@_), 22);
-}
-
-{ # FTP access plus FTP index
- my $t = Tabix->new(-data=>"ftp://ftp.ncbi.nih.gov/1000genomes/ftp/pilot_data/release/2010_03/pilot1/CEU.SRP000031.2010_03.genotypes.vcf.gz",
- -index=>"ftp://ftp.ncbi.nih.gov/1000genomes/ftp/pilot_data/release/2010_03/pilot1/CEU.SRP000031.2010_03.genotypes.vcf.gz.tbi");
- ok($t);
- my $iter = $t->query("19", 10000000, 10100000);
- ok($iter);
- $_ = 0;
- ++$_ while ($t->read($iter));
- is($_, 268);
- @_ = $t->getnames;
- is(scalar(@_), 22);
-}
diff --git a/perl/typemap b/perl/typemap
deleted file mode 100644
index a312f99..0000000
--- a/perl/typemap
+++ /dev/null
@@ -1,3 +0,0 @@
-TYPEMAP
-tabix_t* T_PTROBJ
-ti_iter_t T_PTROBJ
\ No newline at end of file
diff --git a/python/setup.py b/python/setup.py
deleted file mode 100644
index 2771eda..0000000
--- a/python/setup.py
+++ /dev/null
@@ -1,55 +0,0 @@
-#!/usr/bin/env python
-#
-# The MIT License
-#
-# Copyright (c) 2011 Seoul National University.
-#
-# 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: Hyeshik Chang <hyeshik at snu.ac.kr>
-
-from distutils.core import setup, Extension
-
-# Change this to True when you need the knetfile support.
-USE_KNETFILE = False
-
-TABIX_SOURCE_FILES = [
- '../bgzf.c', '../bgzip.c', '../index.c', '../knetfile.c', '../kstring.c'
-]
-
-define_options = [('_FILE_OFFSET_BITS', 64)]
-if USE_KNETFILE:
- define_options.append(('_USE_KNETFILE', 1))
-
-ext_modules = [Extension("tabix", ["tabixmodule.c"] + TABIX_SOURCE_FILES,
- include_dirs=['..'],
- libraries=['z'],
- define_macros=define_options)]
-
-setup (name = 'tabix',
- version = '1.0',
- description = 'Python interface to tabix, a generic indexer '
- 'for TAB-delimited genome position files',
- author = 'Hyeshik Chang',
- author_email = 'hyeshik at snu.ac.kr',
- license = 'MIT',
- ext_modules = ext_modules
-)
diff --git a/python/tabixmodule.c b/python/tabixmodule.c
deleted file mode 100644
index d04d097..0000000
--- a/python/tabixmodule.c
+++ /dev/null
@@ -1,408 +0,0 @@
-/*-
- * The MIT License
- *
- * Copyright (c) 2011 Seoul National University.
- *
- * 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: Hyeshik Chang <hyeshik at snu.ac.kr>
- */
-
-#define PY_SSIZE_T_CLEAN
-#include "Python.h"
-#include "tabix.h"
-
-static PyObject *TabixError;
-
-typedef struct {
- PyObject_HEAD
- tabix_t *tb;
- char *fn;
-} TabixObject;
-
-typedef struct {
- PyObject_HEAD
- TabixObject *tbobj;
- ti_iter_t iter;
-} TabixIteratorObject;
-
-static PyTypeObject Tabix_Type, TabixIterator_Type;
-
-/* --- TabixIterator --------------------------------------------------- */
-
-static PyObject *
-tabixiter_create(TabixObject *parentidx, ti_iter_t iter)
-{
- TabixIteratorObject *self;
-
- self = PyObject_New(TabixIteratorObject, &TabixIterator_Type);
- if (self == NULL)
- return NULL;
-
- Py_INCREF(parentidx);
- self->tbobj = parentidx;
- self->iter = iter;
-
- return (PyObject *)self;
-}
-
-static void
-tabixiter_dealloc(TabixIteratorObject *self)
-{
- ti_iter_destroy(self->iter);
- Py_DECREF(self->tbobj);
- PyObject_Del(self);
-}
-
-static PyObject *
-tabixiter_iter(PyObject *self)
-{
- Py_INCREF(self);
- return self;
-}
-
-#if PY_MAJOR_VERSION < 3
-# define PYOBJECT_FROM_STRING_AND_SIZE PyString_FromStringAndSize
-#else
-# define PYOBJECT_FROM_STRING_AND_SIZE PyUnicode_FromStringAndSize
-#endif
-
-static PyObject *
-tabixiter_iternext(TabixIteratorObject *self)
-{
- const char *chunk;
- int len, i;
-
- chunk = ti_read(self->tbobj->tb, self->iter, &len);
- if (chunk != NULL) {
- PyObject *ret, *column;
- Py_ssize_t colidx;
- const char *ptr, *begin;
-
- ret = PyList_New(0);
- if (ret == NULL)
- return NULL;
-
- colidx = 0;
- ptr = begin = chunk;
- for (i = len; i > 0; i--, ptr++)
- if (*ptr == '\t') {
- column = PYOBJECT_FROM_STRING_AND_SIZE(begin,
- (Py_ssize_t)(ptr - begin));
- if (column == NULL || PyList_Append(ret, column) == -1) {
- Py_DECREF(ret);
- return NULL;
- }
-
- Py_DECREF(column);
- begin = ptr + 1;
- colidx++;
- }
-
- column = PYOBJECT_FROM_STRING_AND_SIZE(begin, (Py_ssize_t)(ptr - begin));
- if (column == NULL || PyList_Append(ret, column) == -1) {
- Py_DECREF(ret);
- return NULL;
- }
- Py_DECREF(column);
-
- return ret;
- }
- else
- return NULL;
-}
-
-static PyMethodDef tabixiter_methods[] = {
- {NULL, NULL} /* sentinel */
-};
-
-static PyTypeObject TabixIterator_Type = {
- PyVarObject_HEAD_INIT(NULL, 0)
- "tabix.TabixIterator", /*tp_name*/
- sizeof(TabixIteratorObject), /*tp_basicsize*/
- 0, /*tp_itemsize*/
- /* methods */
- (destructor)tabixiter_dealloc, /*tp_dealloc*/
- 0, /*tp_print*/
- 0, /*tp_getattr*/
- 0, /*tp_setattr*/
- 0, /*tp_compare*/
- 0, /*tp_repr*/
- 0, /*tp_as_number*/
- 0, /*tp_as_sequence*/
- 0, /*tp_as_mapping*/
- 0, /*tp_hash*/
- 0, /*tp_call*/
- 0, /*tp_str*/
- 0, /*tp_getattro*/
- 0, /*tp_setattro*/
- 0, /*tp_as_buffer*/
- Py_TPFLAGS_DEFAULT, /*tp_flags*/
- 0, /*tp_doc*/
- 0, /*tp_traverse*/
- 0, /*tp_clear*/
- 0, /*tp_richcompare*/
- 0, /*tp_weaklistoffset*/
- tabixiter_iter, /*tp_iter*/
- (iternextfunc)tabixiter_iternext, /*tp_iternext*/
- tabixiter_methods, /*tp_methods*/
- 0, /*tp_members*/
- 0, /*tp_getset*/
- 0, /*tp_base*/
- 0, /*tp_dict*/
- 0, /*tp_descr_get*/
- 0, /*tp_descr_set*/
- 0, /*tp_dictoffset*/
- 0, /*tp_init*/
- 0, /*tp_alloc*/
- 0, /*tp_new*/
- 0, /*tp_free*/
- 0, /*tp_is_gc*/
-};
-
-
-/* --- Tabix ----------------------------------------------------------- */
-
-static PyObject *
-tabix_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
-{
- TabixObject *self;
- const char *fn, *fnidx=NULL;
- static char *kwnames[]={"fn", "fnidx", NULL};
- tabix_t *tb;
-
- if (!PyArg_ParseTupleAndKeywords(args, kwds, "s|z:Tabix",
- kwnames, &fn, &fnidx))
- return NULL;
-
- tb = ti_open(fn, fnidx);
- if (tb == NULL) {
- PyErr_SetString(TabixError, "Can't open the index file.");
- return NULL;
- }
-
- self = (TabixObject *)type->tp_alloc(type, 0);
- if (self == NULL)
- return NULL;
-
- self->tb = tb;
- self->fn = strdup(fn);
-
- return (PyObject *)self;
-}
-
-static void
-tabix_dealloc(TabixObject *self)
-{
- free(self->fn);
- ti_close(self->tb);
- PyObject_Del(self);
-}
-
-static PyObject *
-tabix_query(TabixObject *self, PyObject *args)
-{
- char *name;
- int begin, end;
- ti_iter_t result;
-
- if (!PyArg_ParseTuple(args, "sii:query", &name, &begin, &end))
- return NULL;
-
- result = ti_query(self->tb, name, begin, end);
- if (result == NULL) {
- PyErr_SetString(TabixError, "query failed");
- return NULL;
- }
-
- return tabixiter_create(self, result);
-}
-
-static PyObject *
-tabix_queryi(TabixObject *self, PyObject *args)
-{
- int tid, begin, end;
- ti_iter_t result;
-
- if (!PyArg_ParseTuple(args, "iii:queryi", &tid, &begin, &end))
- return NULL;
-
- result = ti_queryi(self->tb, tid, begin, end);
- if (result == NULL) {
- PyErr_SetString(TabixError, "query failed");
- return NULL;
- }
-
- return tabixiter_create(self, result);
-}
-
-static PyObject *
-tabix_querys(TabixObject *self, PyObject *args)
-{
- const char *reg;
- ti_iter_t result;
-
- if (!PyArg_ParseTuple(args, "s:querys", ®))
- return NULL;
-
- result = ti_querys(self->tb, reg);
- if (result == NULL) {
- PyErr_SetString(TabixError, "query failed");
- return NULL;
- }
-
- return tabixiter_create(self, result);
-}
-
-static PyObject *
-tabix_repr(TabixObject *self)
-{
-#if PY_MAJOR_VERSION < 3
- return PyString_FromFormat("<Tabix fn=\"%s\">", self->fn);
-#else
- return PyUnicode_FromFormat("<Tabix fn=\"%s\">", self->fn);
-#endif
-}
-
-static PyMethodDef tabix_methods[] = {
- {"query", (PyCFunction)tabix_query, METH_VARARGS,
- PyDoc_STR("T.query(name, begin, end) -> iterator")},
- {"queryi", (PyCFunction)tabix_queryi, METH_VARARGS,
- PyDoc_STR("T.queryi(tid, begin, id) -> iterator")},
- {"querys", (PyCFunction)tabix_querys, METH_VARARGS,
- PyDoc_STR("T.querys(region) -> iterator")},
- {NULL, NULL} /* sentinel */
-};
-
-static PyTypeObject Tabix_Type = {
- /* The ob_type field must be initialized in the module init function
- * to be portable to Windows without using C++. */
- PyVarObject_HEAD_INIT(NULL, 0)
- "tabix.Tabix", /*tp_name*/
- sizeof(TabixObject), /*tp_basicsize*/
- 0, /*tp_itemsize*/
- /* methods */
- (destructor)tabix_dealloc, /*tp_dealloc*/
- 0, /*tp_print*/
- 0, /*tp_getattr*/
- 0, /*tp_setattr*/
- 0, /*tp_compare*/
- (reprfunc)tabix_repr, /*tp_repr*/
- 0, /*tp_as_number*/
- 0, /*tp_as_sequence*/
- 0, /*tp_as_mapping*/
- 0, /*tp_hash*/
- 0, /*tp_call*/
- 0, /*tp_str*/
- 0, /*tp_getattro*/
- 0, /*tp_setattro*/
- 0, /*tp_as_buffer*/
- Py_TPFLAGS_DEFAULT, /*tp_flags*/
- 0, /*tp_doc*/
- 0, /*tp_traverse*/
- 0, /*tp_clear*/
- 0, /*tp_richcompare*/
- 0, /*tp_weaklistoffset*/
- 0, /*tp_iter*/
- 0, /*tp_iternext*/
- tabix_methods, /*tp_methods*/
- 0, /*tp_members*/
- 0, /*tp_getset*/
- 0, /*tp_base*/
- 0, /*tp_dict*/
- 0, /*tp_descr_get*/
- 0, /*tp_descr_set*/
- 0, /*tp_dictoffset*/
- 0, /*tp_init*/
- 0, /*tp_alloc*/
- (newfunc)tabix_new, /*tp_new*/
- 0, /*tp_free*/
- 0, /*tp_is_gc*/
-};
-/* --------------------------------------------------------------------- */
-
-static PyMethodDef tabix_functions[] = {
- {NULL, NULL} /* sentinel */
-};
-
-PyDoc_STRVAR(module_doc,
-"Python interface to tabix, Heng Li's generic indexer for TAB-delimited "
-"genome position filesThis is a template module just for instruction.");
-
-#if PY_MAJOR_VERSION >= 3
-static struct PyModuleDef tabixmodule = {
- PyModuleDef_HEAD_INIT,
- "tabix",
- module_doc,
- -1,
- tabix_functions,
- NULL,
- NULL,
- NULL,
- NULL
-};
-#endif
-
-#if PY_MAJOR_VERSION < 3
-PyMODINIT_FUNC inittabix(void)
-#else
-PyMODINIT_FUNC PyInit_tabix(void)
-#endif
-{
- PyObject *m;
-
- if (PyType_Ready(&Tabix_Type) < 0)
- goto fail;
- if (PyType_Ready(&TabixIterator_Type) < 0)
- goto fail;
-
-#if PY_MAJOR_VERSION < 3
- m = Py_InitModule3("tabix", tabix_functions, module_doc);
-#else
- m = PyModule_Create(&tabixmodule);
-#endif
- if (m == NULL)
- goto fail;
-
- if (TabixError == NULL) {
- TabixError = PyErr_NewException("tabix.error", NULL, NULL);
- if (TabixError == NULL)
- goto fail;
- }
- Py_INCREF(TabixError);
- PyModule_AddObject(m, "error", TabixError);
-
- PyModule_AddObject(m, "Tabix", (PyObject *)&Tabix_Type);
- PyModule_AddObject(m, "TabixIterator", (PyObject *)&TabixIterator_Type);
-
-#if PY_MAJOR_VERSION >= 3
- return m;
-#endif
-
- fail:
-#if PY_MAJOR_VERSION < 3
- return;
-#else
- return NULL;
-#endif
-}
diff --git a/python/test.py b/python/test.py
deleted file mode 100755
index 804a565..0000000
--- a/python/test.py
+++ /dev/null
@@ -1,91 +0,0 @@
-#
-# The MIT License
-#
-# Copyright (c) 2011 Seoul National University.
-#
-# 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: Hyeshik Chang <hyeshik at snu.ac.kr>
-
-import unittest
-import random
-import gzip
-import tabix
-
-EXAMPLEFILE = '../example.gtf.gz'
-
-def load_example_regions(path):
- alldata = []
- for line in gzip.GzipFile(EXAMPLEFILE):
- fields = line.decode('ascii')[:-1].split('\t')
- seqid = fields[0]
- begin = int(fields[3])
- end = int(fields[4])
- alldata.append((seqid, begin, end, fields[:7]))
-
- return alldata
-
-def does_overlap(A, B, C, D):
- return (A <= D <= B) or (C <= B <= D)
-
-def sample_test_dataset(regions, ntests):
- seqids = [seqid for seqid, _, _, _ in regions]
- lowerbound = max(0, min(begin for _, begin, _, _ in regions) - 1000)
- upperbound = max(end for _, _, end, _ in regions) + 1000
-
- tests = []
- for i in range(ntests):
- seqid = random.choice(seqids)
- low = random.randrange(lowerbound, upperbound)
- high = random.randrange(low, upperbound)
-
- # for 1-based both-end inclusive intervals
- matches = [info for seq, begin, end, info in regions
- if seqid == seq and does_overlap(begin, end, low, high)]
-
- tests.append((seqid, low, high, matches))
-
- return tests
-
-def tbresult2excerpt(tbmatches):
- return [fields[:7] for fields in tbmatches]
-
-class TabixTest(unittest.TestCase):
- regions = load_example_regions(EXAMPLEFILE)
- testset = sample_test_dataset(regions, 500)
-
- def setUp(self):
- self.tb = tabix.Tabix(EXAMPLEFILE)
-
- def testQuery(self):
- for seqid, low, high, matches in self.testset:
- tbresult = tbresult2excerpt(self.tb.query(seqid, low, high))
- self.assertEqual(tbresult, matches)
-
- def testQueryS(self):
- for seqid, low, high, matches in self.testset:
- tbresult = tbresult2excerpt(self.tb.querys('%s:%d-%d' %
- (seqid, low, high)))
- self.assertEqual(tbresult, matches)
-
-
-if __name__ == '__main__':
- unittest.main()
diff --git a/tabix.1 b/tabix.1
deleted file mode 100644
index 1bd9533..0000000
--- a/tabix.1
+++ /dev/null
@@ -1,132 +0,0 @@
-.TH tabix 1 "11 May 2010" "tabix-0.2.0" "Bioinformatics tools"
-.SH NAME
-.PP
-bgzip - Block compression/decompression utility
-.PP
-tabix - Generic indexer for TAB-delimited genome position files
-.SH SYNOPSIS
-.PP
-.B bgzip
-.RB [ \-cdhB ]
-.RB [ \-b
-.IR virtualOffset ]
-.RB [ \-s
-.IR size ]
-.RI [ file ]
-.PP
-.B tabix
-.RB [ \-0lf ]
-.RB [ \-p
-.R gff|bed|sam|vcf]
-.RB [ \-s
-.IR seqCol ]
-.RB [ \-b
-.IR begCol ]
-.RB [ \-e
-.IR endCol ]
-.RB [ \-S
-.IR lineSkip ]
-.RB [ \-c
-.IR metaChar ]
-.I in.tab.bgz
-.RI [ "region1 " [ "region2 " [ ... "]]]"
-
-.SH DESCRIPTION
-.PP
-Tabix indexes a TAB-delimited genome position file
-.I in.tab.bgz
-and creates an index file
-.I in.tab.bgz.tbi
-when
-.I region
-is absent from the command-line. The input data file must be position
-sorted and compressed by
-.B bgzip
-which has a
-.BR gzip (1)
-like interface. After indexing, tabix is able to quickly retrieve data
-lines overlapping
-.I regions
-specified in the format "chr:beginPos-endPos". Fast data retrieval also
-works over network if URI is given as a file name and in this case the
-index file will be downloaded if it is not present locally.
-
-.SH OPTIONS OF TABIX
-.TP 10
-.BI "-p " STR
-Input format for indexing. Valid values are: gff, bed, sam, vcf and
-psltab. This option should not be applied together with any of
-.BR \-s ", " \-b ", " \-e ", " \-c " and " \-0 ;
-it is not used for data retrieval because this setting is stored in
-the index file. [gff]
-.TP
-.BI "-s " INT
-Column of sequence name. Option
-.BR \-s ", " \-b ", " \-e ", " \-S ", " \-c " and " \-0
-are all stored in the index file and thus not used in data retrieval. [1]
-.TP
-.BI "-b " INT
-Column of start chromosomal position. [4]
-.TP
-.BI "-e " INT
-Column of end chromosomal position. The end column can be the same as the
-start column. [5]
-.TP
-.BI "-S " INT
-Skip first INT lines in the data file. [0]
-.TP
-.BI "-c " CHAR
-Skip lines started with character CHAR. [#]
-.TP
-.B -0
-Specify that the position in the data file is 0-based (e.g. UCSC files)
-rather than 1-based.
-.TP
-.B -h
-Print the header/meta lines.
-.TP
-.B -B
-The second argument is a BED file. When this option is in use, the input
-file may not be sorted or indexed. The entire input will be read sequentially. Nonetheless,
-with this option, the format of the input must be specificed correctly on the command line.
-.TP
-.B -f
-Force to overwrite the index file if it is present.
-.TP
-.B -l
-List the sequence names stored in the index file.
-.RE
-
-.SH EXAMPLE
-(grep ^"#" in.gff; grep -v ^"#" in.gff | sort -k1,1 -k4,4n) | bgzip > sorted.gff.gz;
-
-tabix -p gff sorted.gff.gz;
-
-tabix sorted.gff.gz chr1:10,000,000-20,000,000;
-
-.SH NOTES
-It is straightforward to achieve overlap queries using the standard
-B-tree index (with or without binning) implemented in all SQL databases,
-or the R-tree index in PostgreSQL and Oracle. But there are still many
-reasons to use tabix. Firstly, tabix directly works with a lot of widely
-used TAB-delimited formats such as GFF/GTF and BED. We do not need to
-design database schema or specialized binary formats. Data do not need
-to be duplicated in different formats, either. Secondly, tabix works on
-compressed data files while most SQL databases do not. The GenCode
-annotation GTF can be compressed down to 4%. Thirdly, tabix is
-fast. The same indexing algorithm is known to work efficiently for an
-alignment with a few billion short reads. SQL databases probably cannot
-easily handle data at this scale. Last but not the least, tabix supports
-remote data retrieval. One can put the data file and the index at an FTP
-or HTTP server, and other users or even web services will be able to get
-a slice without downloading the entire file.
-
-.SH AUTHOR
-.PP
-Tabix was written by Heng Li. The BGZF library was originally
-implemented by Bob Handsaker and modified by Heng Li for remote file
-access and in-memory caching.
-
-.SH SEE ALSO
-.PP
-.BR samtools (1)
diff --git a/tabix.cpp b/tabix.cpp
index 1d75996..73ff3e1 100644
--- a/tabix.cpp
+++ b/tabix.cpp
@@ -3,12 +3,13 @@
Tabix::Tabix(void) { }
Tabix::Tabix(string& file) {
+ has_jumped = false;
filename = file;
const char* cfilename = file.c_str();
struct stat stat_tbi,stat_vcf;
char *fnidx = (char*) calloc(strlen(cfilename) + 5, 1);
strcat(strcpy(fnidx, cfilename), ".tbi");
- if ( bgzf_check_bgzf(cfilename)!=1 )
+ if ( bgzf_is_bgzf(cfilename)!=1 )
{
cerr << "[tabix++] was bgzip used to compress this file? " << file << endl;
free(fnidx);
@@ -25,65 +26,85 @@ Tabix::Tabix(string& file) {
}
free(fnidx);
- if ((t = ti_open(cfilename, 0)) == 0) {
+ if ((fn = hts_open(cfilename, "r")) == 0) {
cerr << "[tabix++] fail to open the data file." << endl;
exit(1);
}
- if (ti_lazy_index_load(t) < 0) {
+ if ((tbx = tbx_index_load(cfilename)) == NULL) {
cerr << "[tabix++] failed to load the index file." << endl;
exit(1);
}
- idxconf = ti_get_conf(t->idx);
+ int nseq;
+ const char** seq = tbx_seqnames(tbx, &nseq);
+ for (int i=0; i<nseq; i++) {
+ chroms.push_back(seq[i]);
+ }
+ free(seq);
+
+ idxconf = &tbx_conf_vcf;
// set up the iterator, defaults to the beginning
- iter = ti_query(t, 0, 0, 0);
+ current_chrom = chroms.begin();
+ iter = tbx_itr_querys(tbx, current_chrom->c_str());
}
Tabix::~Tabix(void) {
- ti_iter_destroy(iter);
- ti_close(t);
+ tbx_itr_destroy(iter);
+ tbx_destroy(tbx);
}
-
void Tabix::getHeader(string& header) {
header.clear();
- ti_iter_destroy(iter);
- iter = ti_query(t, 0, 0, 0);
- const char* s;
- int len;
- while ((s = ti_read(t, iter, &len)) != 0) {
- if ((int)(*s) != idxconf->meta_char) {
- firstline = string(s); // stash this line
+ kstring_t str = {0,0,0};
+ while ( hts_getline(fn, KS_SEP_LINE, &str) >= 0 ) {
+ if ( !str.l || str.s[0]!=tbx->conf.meta_char ) {
break;
} else {
- header += string(s);
+ header += string(str.s);
header += "\n";
}
}
+ // set back to start
+ current_chrom = chroms.begin();
+ if (iter) tbx_itr_destroy(iter);
+ iter = tbx_itr_querys(tbx, current_chrom->c_str());
}
bool Tabix::setRegion(string& region) {
- if (ti_parse_region(t->idx, region.c_str(), &tid, &beg, &end) == 0) {
- firstline.clear();
- ti_iter_destroy(iter);
- iter = ti_queryi(t, tid, beg, end);
- return true;
- } else return false;
+ tbx_itr_destroy(iter);
+ iter = tbx_itr_querys(tbx, region.c_str());
+ has_jumped = true;
+ return true;
}
bool Tabix::getNextLine(string& line) {
- const char* s;
- int len;
- if (!firstline.empty()) {
- line = firstline; // recovers line read if header is parsed
- firstline.clear();
- return true;
+ kstring_t str = {0,0,0};
+ if (has_jumped) {
+ if (iter && tbx_itr_next(fn, tbx, iter, &str) >= 0) {
+ line = string(str.s);
+ return true;
+ } else return false;
+ } else { // step through all sequences in the file
+ // we've never jumped, so read everything
+ if (iter && tbx_itr_next(fn, tbx, iter, &str) >= 0) {
+ line = string(str.s);
+ return true;
+ } else {
+ ++current_chrom;
+ while (current_chrom != chroms.end()) {
+ tbx_itr_destroy(iter);
+ iter = tbx_itr_querys(tbx, current_chrom->c_str());
+ if (iter && tbx_itr_next(fn, tbx, iter, &str) >= 0) {
+ line = string(str.s);
+ return true;
+ } else {
+ ++current_chrom;
+ }
+ }
+ return false;
+ }
}
- if ((s = ti_read(t, iter, &len)) != 0) {
- line = string(s);
- return true;
- } else return false;
}
diff --git a/tabix.h b/tabix.h
deleted file mode 100644
index 7b4497a..0000000
--- a/tabix.h
+++ /dev/null
@@ -1,145 +0,0 @@
-/* The MIT License
-
- Copyright (c) 2009 Genome Research Ltd (GRL), 2010 Broad Institute
-
- 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 at live.co.uk> */
-
-#ifndef __TABIDX_H
-#define __TABIDX_H
-
-#include <stdint.h>
-#include "kstring.h"
-#include "bgzf.h"
-
-#define TI_PRESET_GENERIC 0
-#define TI_PRESET_SAM 1
-#define TI_PRESET_VCF 2
-
-#define TI_FLAG_UCSC 0x10000
-
-typedef int (*ti_fetch_f)(int l, const char *s, void *data);
-
-struct __ti_index_t;
-typedef struct __ti_index_t ti_index_t;
-
-struct __ti_iter_t;
-typedef struct __ti_iter_t *ti_iter_t;
-
-typedef struct {
- BGZF *fp;
- ti_index_t *idx;
- char *fn, *fnidx;
-} tabix_t;
-
-typedef struct {
- int32_t preset;
- int32_t sc, bc, ec; // seq col., beg col. and end col.
- int32_t meta_char, line_skip;
-} ti_conf_t;
-
-typedef struct {
- int beg, end;
- char *ss, *se;
-} ti_interval_t;
-
-extern ti_conf_t ti_conf_gff, ti_conf_bed, ti_conf_psltbl, ti_conf_vcf, ti_conf_sam; // preset
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
- /*******************
- * High-level APIs *
- *******************/
-
- tabix_t *ti_open(const char *fn, const char *fnidx);
- int ti_lazy_index_load(tabix_t *t);
- void ti_close(tabix_t *t);
- ti_iter_t ti_query(tabix_t *t, const char *name, int beg, int end);
- ti_iter_t ti_queryi(tabix_t *t, int tid, int beg, int end);
- ti_iter_t ti_querys(tabix_t *t, const char *reg);
- const char *ti_read(tabix_t *t, ti_iter_t iter, int *len);
-
- /* Destroy the iterator */
- void ti_iter_destroy(ti_iter_t iter);
-
- /* Get the list of sequence names. Each "char*" pointer points to a
- * internal member of the index, so DO NOT modify the returned
- * pointer; otherwise the index will be corrupted. The returned
- * pointer should be freed by a single free() call by the routine
- * calling this function. The number of sequences is returned at *n. */
- const char **ti_seqname(const ti_index_t *idx, int *n);
-
- /******************
- * Low-level APIs *
- ******************/
-
- /* Build the index for file <fn>. File <fn>.tbi will be generated
- * and overwrite the file of the same name. Return -1 on failure. */
- int ti_index_build(const char *fn, const ti_conf_t *conf);
-
- /* Load the index from file <fn>.tbi. If <fn> is a URL and the index
- * file is not in the working directory, <fn>.tbi will be
- * downloaded. Return NULL on failure. */
- ti_index_t *ti_index_load(const char *fn);
-
- ti_index_t *ti_index_load_local(const char *fnidx);
-
- /* Destroy the index */
- void ti_index_destroy(ti_index_t *idx);
-
- /* Parse a region like: chr2, chr2:100, chr2:100-200. Return -1 on failure. */
- int ti_parse_region(const ti_index_t *idx, const char *str, int *tid, int *begin, int *end);
-
- int ti_get_tid(const ti_index_t *idx, const char *name);
-
- /* Get the iterator pointing to the first record at the current file
- * position. If the file is just openned, the iterator points to the
- * first record in the file. */
- ti_iter_t ti_iter_first(void);
-
- /* Get the iterator pointing to the first record in region tid:beg-end */
- ti_iter_t ti_iter_query(const ti_index_t *idx, int tid, int beg, int end);
-
- /* Get the data line pointed by the iterator and iterate to the next record. */
- const char *ti_iter_read(BGZF *fp, ti_iter_t iter, int *len);
-
- const ti_conf_t *ti_get_conf(ti_index_t *idx);
- int ti_get_intv(const ti_conf_t *conf, int len, char *line, ti_interval_t *intv);
-
- /*******************
- * Deprecated APIs *
- *******************/
-
- /* The callback version for random access */
- int ti_fetch(BGZF *fp, const ti_index_t *idx, int tid, int beg, int end, void *data, ti_fetch_f func);
-
- /* Read one line. */
- int ti_readline(BGZF *fp, kstring_t *str);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
diff --git a/tabix.hpp b/tabix.hpp
index 7b67a6f..e9ab48d 100644
--- a/tabix.hpp
+++ b/tabix.hpp
@@ -1,24 +1,31 @@
#include <string>
#include <stdlib.h>
#include <sys/stat.h>
-#include "bgzf.h"
-#include "tabix.h"
+#include "htslib/bgzf.h"
+#include "htslib/tbx.h"
+#include "htslib/kseq.h"
#include <iostream>
+#include <cstring>
+#include <vector>
using namespace std;
class Tabix {
- tabix_t *t;
- ti_iter_t iter;
- const ti_conf_t *idxconf;
+ htsFile* fn;
+ tbx_t* tbx;
+ hts_itr_t* iter;
+ const tbx_conf_t *idxconf;
int tid, beg, end;
string firstline;
+ bool has_jumped;
+ vector<string>::iterator current_chrom;
public:
string filename;
+ vector<string> chroms;
Tabix(void);
Tabix(string& file);
diff --git a/tabix.py b/tabix.py
deleted file mode 100755
index 8b5c7c8..0000000
--- a/tabix.py
+++ /dev/null
@@ -1,87 +0,0 @@
-#!/usr/bin/env python
-
-# Author: Heng Li and Aaron Quinlan
-# License: MIT/X11
-
-import sys
-from ctypes import *
-from ctypes.util import find_library
-import glob, platform
-
-def load_shared_library(lib, _path='.', ver='*'):
- """Search for and load the tabix library. The
- expectation is that the library is located in
- the current directory (ie. "./")
- """
- # find from the system path
- path = find_library(lib)
- if (path == None): # if fail, search in the custom directory
- s = platform.system()
- if (s == 'Darwin'): suf = ver+'.dylib'
- elif (s == 'Linux'): suf = '.so'+ver
- candidates = glob.glob(_path+'/lib'+lib+suf);
- if (len(candidates) == 1): path = candidates[0]
- else: return None
- cdll.LoadLibrary(path)
- return CDLL(path)
-
-def tabix_init():
- """Initialize and return a tabix reader object
- for subsequent tabix_get() calls.
- """
- tabix = load_shared_library('tabix')
- if (tabix == None): return None
- tabix.ti_read.restype = c_char_p
- # on Mac OS X 10.6, the following declarations are required.
- tabix.ti_open.restype = c_void_p
- tabix.ti_querys.argtypes = [c_void_p, c_char_p]
- tabix.ti_querys.restype = c_void_p
- tabix.ti_query.argtypes = [c_void_p, c_char_p, c_int, c_int]
- tabix.ti_query.restype = c_void_p
- tabix.ti_read.argtypes = [c_void_p, c_void_p, c_void_p]
- tabix.ti_iter_destroy.argtypes = [c_void_p]
- tabix.ti_close.argtypes = [c_void_p]
- # FIXME: explicit declarations for APIs not used in this script
- return tabix
-
-# OOP interface
-class Tabix:
- def __init__(self, fn, fnidx=0):
- self.tabix = tabix_init();
- if (self.tabix == None):
- sys.stderr.write("[Tabix] Please make sure the shared library is compiled and available.\n")
- return
- self.fp = self.tabix.ti_open(fn, fnidx);
-
- def __del__(self):
- if (self.tabix): self.tabix.ti_close(self.fp)
-
- def fetch(self, chr, start=-1, end=-1):
- """Generator function that will yield each interval
- within the requested range from the requested file.
- """
- if (self.tabix == None): return
- if (start < 0): iter = self.tabix.ti_querys(self.fp, chr) # chr looks like: "chr2:1,000-2,000" or "chr2"
- else: iter = self.tabix.ti_query(self.fp, chr, start, end) # chr must be a sequence name
- if (iter == None):
- sys.stderr.write("[Tabix] Malformatted query or wrong sequence name.\n")
- return
- while (1): # iterate
- s = self.tabix.ti_read(self.fp, iter, 0)
- if (s == None): break
- yield s
- self.tabix.ti_iter_destroy(iter)
-
-# command-line interface
-def main():
- if (len(sys.argv) < 3):
- sys.stderr.write("Usage: tabix.py <in.gz> <reg>\n")
- sys.exit(1)
-
- # report the features in the requested interval
- tabix = Tabix(sys.argv[1])
- for line in tabix.fetch(sys.argv[2]):
- print line
-
-if __name__ == '__main__':
- main()
diff --git a/tabix.tex b/tabix.tex
deleted file mode 100644
index 669d673..0000000
--- a/tabix.tex
+++ /dev/null
@@ -1,121 +0,0 @@
-\documentclass[10pt]{article}
-\usepackage{color}
-\definecolor{gray}{rgb}{0.7,0.7,0.7}
-
-\setlength{\topmargin}{0.0cm}
-\setlength{\textheight}{21.5cm}
-\setlength{\oddsidemargin}{0cm}
-\setlength{\textwidth}{16.5cm}
-\setlength{\columnsep}{0.6cm}
-
-\begin{document}
-
-\title{The Tabix index file format}
-\author{Heng Li}
-\date{}
-
-\maketitle
-
-\begin{center}
-\begin{tabular}{|l|l|l|l|l|l|l|}
-\hline
-\multicolumn{4}{|c|}{\bf Field} & \multicolumn{1}{c|}{\bf Descrption} & \multicolumn{1}{c|}{\bf Type} & \multicolumn{1}{c|}{\bf Value} \\
-\hline\hline
-\multicolumn{4}{|l|}{\tt magic} & Magic string & {\tt char[4]} & TBI$\backslash$1 \\
-\hline
-\multicolumn{4}{|l|}{\tt n\_ref} & \# sequences & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt format} & Format (0: generic; 1: SAM; 2: VCF) & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt col\_seq} & Column for the sequence name & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt col\_beg} & Column for the start of a region & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt col\_end} & Column for the end of a region & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt meta} & Leading character for comment lines & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt skip} & \# lines to skip at the beginning & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt l\_nm} & Length of concatenated sequence names & {\tt int32\_t} & \\
-\hline
-\multicolumn{4}{|l|}{\tt names} & Concatenated names, each zero terminated & {\tt char[l\_nm]} & \\
-\hline
-\multicolumn{7}{|c|}{\textcolor{gray}{\it List of indices (n=n\_ref)}}\\
-\cline{2-7}
-\hspace{0.1cm} & \multicolumn{3}{l|}{\tt n\_bin} & \# distinct bins (for the binning index) & {\tt int32\_t} & \\
-\cline{2-7}
- & \multicolumn{6}{c|}{\textcolor{gray}{\it List of distinct bins (n=n\_bin)}} \\
-\cline{3-7}
- & \hspace{0.1cm} & \multicolumn{2}{l|}{\tt bin} & Distinct bin number & {\tt uint32\_t} & \\
-\cline{3-7}
- & & \multicolumn{2}{l|}{\tt n\_chunk} & \# chunks & {\tt int32\_t} & \\
-\cline{3-7}
- & & \multicolumn{5}{c|}{\textcolor{gray}{\it List of chunks (n=n\_chunk)}} \\
-\cline{4-7}
- & & \hspace{0.1cm} & {\tt cnk\_beg} & Virtual file offset of the start of the chunk & {\tt uint64\_t} & \\
-\cline{4-7}
- & & & {\tt cnk\_end} & Virtual file offset of the end of the chunk & {\tt uint64\_t} & \\
-\cline{2-7}
- & \multicolumn{3}{l|}{\tt n\_intv} & \# 16kb intervals (for the linear index) & {\tt int32\_t} & \\
-\cline{2-7}
- & \multicolumn{6}{c|}{\textcolor{gray}{\it List of distinct intervals (n=n\_intv)}} \\
-\cline{3-7}
- & & \multicolumn{2}{l|}{\tt ioff} & File offset of the first record in the interval & {\tt uint64\_t} & \\
-\hline
-\end{tabular}
-\end{center}
-
-{\bf Notes:}
-
-\begin{itemize}
-\item The index file is BGZF compressed.
-\item All integers are little-endian.
-\item When {\tt (format\&0x10000)} is true, the coordinate follows the
- {\tt BED} rule (i.e. half-closed-half-open and zero based); otherwise,
- the coordinate follows the {\tt GFF} rule (closed and one based).
-\item For the SAM format, the end of a region equals {\tt POS} plus the
- reference length in the alignment, inferred from {\tt CIGAR}. For the
- VCF format, the end of a region equals {\tt POS} plus the size of the
- deletion.
-\item Field {\tt col\_beg} may equal {\tt col\_end}, and in this case,
- the end of a region is {\tt end}={\tt beg+1}.
-\item Example. For {\tt GFF}, {\tt format}=0, {\tt col\_seq}=1, {\tt
- col\_beg}=4, {\tt col\_end}=5, {\tt meta}=`{\tt \#}' and {\tt
- skip}=0. For {\tt BED}, {\tt format}=0x10000, {\tt col\_seq}=1, {\tt
- col\_beg}=2, {\tt col\_end}=3, {\tt meta}=`{\tt \#}' and {\tt
- skip}=0.
-\item Given a zero-based, half-closed and half-open region {\tt
- [beg,end)}, the {\tt bin} number is calculated with the following C
- function:
-\begin{verbatim}
-int reg2bin(int beg, int end) {
- --end;
- if (beg>>14 == end>>14) return ((1<<15)-1)/7 + (beg>>14);
- if (beg>>17 == end>>17) return ((1<<12)-1)/7 + (beg>>17);
- if (beg>>20 == end>>20) return ((1<<9)-1)/7 + (beg>>20);
- if (beg>>23 == end>>23) return ((1<<6)-1)/7 + (beg>>23);
- if (beg>>26 == end>>26) return ((1<<3)-1)/7 + (beg>>26);
- return 0;
-}
-\end{verbatim}
-\item The list of bins that may overlap a region {\tt [beg,end)} can be
- obtained with the following C function.
-\begin{verbatim}
-#define MAX_BIN (((1<<18)-1)/7)
-int reg2bins(int rbeg, int rend, uint16_t list[MAX_BIN])
-{
- int i = 0, k;
- --rend;
- list[i++] = 0;
- for (k = 1 + (rbeg>>26); k <= 1 + (rend>>26); ++k) list[i++] = k;
- for (k = 9 + (rbeg>>23); k <= 9 + (rend>>23); ++k) list[i++] = k;
- for (k = 73 + (rbeg>>20); k <= 73 + (rend>>20); ++k) list[i++] = k;
- for (k = 585 + (rbeg>>17); k <= 585 + (rend>>17); ++k) list[i++] = k;
- for (k = 4681 + (rbeg>>14); k <= 4681 + (rend>>14); ++k) list[i++] = k;
- return i; // #elements in list[]
-}
-\end{verbatim}
-\end{itemize}
-
-\end{document}
\ No newline at end of file
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/libtabixpp.git
More information about the debian-med-commit
mailing list