diff --git a/Makefile b/Makefile index 1fb9b73..4e35554 100644 --- a/Makefile +++ b/Makefile @@ -1,56 +1,74 @@ CC= gcc -CFLAGS= -g -Wall -O2 -DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -DBGZF_CACHE -LOBJS= bgzf.o kstring.o bam_aux.o bam.o bam_import.o bam_index.o sam_header.o bedutil.o commons.o -PROG= bamdst -INCLUDES= -Isamlib/ -I. -SUBDIRS= . samlib -LIBPATH= -L. +CFLAGS= -g -Wall -O2 -std=gnu99 +DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE + +# Source files +SOURCES= xamdst.c kstring.c bedutil.c commons.c +PROG= xamdst + +# ---- htslib configuration --------------------------------------------------- +# Resolution order: +# 1. HTSLIB_DIR= explicit override (a built htslib source tree) +# 2. system htslib detected via pkg-config +# 3. sibling ./htslib local source tree -- libhts.a is auto-built if missing +# 4. otherwise error +# +# If static linking fails with "cannot find -lbz2/-llzma", your htslib was +# built without them -- override: +# make HTSLIB_DEPS="-lz -lpthread -lm" +HTSLIB_DIR ?= + +ifeq ($(HTSLIB_DIR),) + HAVE_PKG_HTSLIB := $(shell pkg-config --exists htslib 2>/dev/null && echo yes) + ifneq ($(HAVE_PKG_HTSLIB),) + # 2. system htslib (pkg-config) + HTSLIB_CFLAGS := $(shell pkg-config --cflags htslib) + HTSLIB_LIBS := $(shell pkg-config --libs htslib) + else + ifneq ($(wildcard htslib/htslib/hts.h),) + # 3. sibling ./htslib source tree (libhts.a auto-built below) + HTSLIB_CFLAGS := -Ihtslib + HTSLIB_LIBS := htslib/libhts.a + HTSLIB_DEPS ?= -lz -lbz2 -llzma -lpthread -lm + HTSLIB_BUILD := htslib/libhts.a + else + # 4. not found + $(error htslib not found. \ +Install htslib system-wide (then pkg-config detects it), \ +or clone it as a sibling: `git clone https://github.com/samtools/htslib htslib`, \ +or pass HTSLIB_DIR=/path/to/htslib) + endif + endif +else + # 1. explicit override + HTSLIB_CFLAGS := -I$(HTSLIB_DIR) + HTSLIB_LIBS := $(HTSLIB_DIR)/libhts.a + HTSLIB_DEPS ?= -lz -lbz2 -llzma -lpthread -lm + HTSLIB_BUILD := $(HTSLIB_DIR)/libhts.a +endif + +INCLUDES= -I. $(HTSLIB_CFLAGS) +LIBS= $(HTSLIB_LIBS) $(HTSLIB_DEPS) -lz -lpthread -lm .SUFFIXES:.c .o -.PHONY: all +.PHONY: all clean + +# auto-build the local htslib (handles a git clone: autoreconf + configure + make) +ifdef HTSLIB_BUILD +$(HTSLIB_BUILD): + cd $(dir $(HTSLIB_BUILD)) && \ + (test -f configure || autoreconf -i) && \ + (test -f config.status || ./configure) && \ + $(MAKE) +endif .c.o: $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $< -o $@ -all:clean $(PROG) - -.PHONY:all clean - -lib:libbam.a - -libbam.a:$(LOBJS) - $(AR) -csru $@ $(LOBJS) - -bamdst:lib $(AOBJS) samlib/bam.h - $(CC) $(CFLAGS) -o $@ $(AOBJS) $(LDFLAGS) bamdst.c $(LIBPATH) $(INCLUDES) -lm -lbam -lz - -bgzf.o:bgzf.c bgzf.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) bgzf.c -o $@ - -kstring.o:kstring.c kstring.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) kstring.c -o $@ - -bam.o:samlib/bam.h samlib/bam_endian.h kstring.h samlib/sam_header.h samlib/bam.c - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) samlib/bam.c -o $@ - -bam_import.o:samlib/bam.h kseq.h khash.h samlib/bam_import.c - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) samlib/bam_import.c -o $@ - -bam_index.o:samlib/bam.h khash.h ksort.h samlib/bam_endian.h samlib/bam_index.c - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) samlib/bam_index.c -o $@ - -sam_header.o:samlib/sam_header.h khash.h samlib/sam_header.c - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) samlib/sam_header.c -o $@ - -bam_aux.o:samlib/bam.h samlib/bam_aux.c khash.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) samlib/bam_aux.c -o $@ - -commons.o:commons.c commons.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) commons.c -o $@ +all:clean $(PROG) -bedutil.o:bedutil.c bedutil.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) bedutil.c -o $@ +xamdst: $(SOURCES:.c=.o) $(HTSLIB_BUILD) + $(CC) $(CFLAGS) -o $@ $(SOURCES:.c=.o) $(INCLUDES) $(LIBS) clean: rm -fr gmon.out *.o a.out *.exe *.dSYM $(PROG) *~ *.a target.dep *.plot *.report *.tsv.gz uncover.bed diff --git a/README.md b/README.md index 99efe9c..71cb5f2 100644 --- a/README.md +++ b/README.md @@ -1,223 +1,115 @@ -# bamdst -- a BAM Depth Stat Tool +# xamdst -Bamdst is a lightweight tool to stat the depth coverage of target regions of bam file(s). +A fast depth-coverage statistic tool for **BAM/CRAM/SAM** files over target regions. +A hard fork of [bamdst](https://github.com/shiquan/bamdst), refactored on +[HTSlib](https://github.com/samtools/htslib) — adding CRAM support, multi-threaded +I/O, and a JSON report. -Bam file(s) should be properly sorted, and the probe file (bed file) and the output dir +Input files must be **coordinate-sorted**. The probe file and output directory are +mandatory. -must be specified in the first place. +## Features -## Install -``` -git clone https://github.com/shiquan/bamdst -cd bamdst -make -``` - -Or by conda. Thanks @xdgene -``` -conda install xdgene::bamdst -``` - - -## USAGE - -Normal: - - bamdst -p -o ./ in1.bam - -Pipeline mode: - - samtools view in1.bam -u | bamdst -p x.bed -o ./ - - -## PARAMETERS - --o / --outdir [dir] - -set the output dir [mandatory] - --p / --bed [file] - -the probe or captured target region file, these regions will be merged first [mandatory] - -## OPTIONAL PARAMETERS - --f / --flank [num] - -if you want calculate the coverage of flank region, set this value, default is 200 - ---maxdepth [num] - -for some projects, the depths of sepcial region are very high, if you don't want show - -these unnormal depths in cumulation distrbution file, set the cutoff value to filter them. - -default is 0 (no filter). - ---cutoffdepth [num] - -for some projects, people care about the coverage of specified depth, like 10000x etc. - -bamdst just calculate the coverage of 0x, 4x, 10x, 30x, 100x, so you can set this value - -to show the specified coverage in the coverage.report file. Default is 0. - ---isize [num] - -for bad mapped paired reads, the inferred insert size is very huge. So set a cutoff - -value for reasonal visual purpose. Default is 2000. - ---uncover [num] - -set this cutoff value for calculate the bad covered region. Default is <5. - ---use_rmdup (an invalid parament since v1.0.0 ) - -Use rmdup depth instead of cover depth to calculate the coverage of target regions and - -so on. +- BAM / CRAM / SAM input, single file or multiple files merged +- Multi-threaded compression & decompression (HTSlib native) +- Per-position, per-region, per-chromosome and cumulative depth statistics +- `raw`, `rmdup` (dedup + primary + mapQ filter) and `coverage` (CIGAR-aware) depths +- JSON report for programmatic access +- `depth.tsv.gz` and `region.tsv.gz` are bgzipped and tabix-ready -## OUTPUT FILES +## Dependencies -Seven files will be created in the output direction. There are: +- **htslib**, **zlib**, **pthread** --**coverage.report** - --**cumu.plot** - --**insert.plot** - --**chromosome.report** - --**region.tsv.gz** - --**depth.tsv.gz** - --**uncover.bed** - -## DETAILS of each file - -**coverage.report** - -This file contains all the coverage information of target and - -flank region, and reads stat information of the input file(s). - -Here is the full details of each entry. - - [Total] Raw Reads (All reads) // All reads in the bam file(s). - [Total] QC Fail reads // Reads number failed QC, this flag is marked by other software,like bwa. See flag in the bam structure. - [Total] Raw Data(Mb) // Total reads data in the bam file(s). - [Total] Paired Reads // Paired reads numbers. - [Total] Mapped Reads // Mapped reads numbers. - [Total] Fraction of Mapped Reads // Ratio of mapped reads against raw reads. - [Total] Mapped Data(Mb) // Mapped data in the bam file(s). - [Total] Fraction of Mapped Data(Mb) // Ratio of mapped data against raw data. - [Total] Properly paired // Paired reads with properly insert size. See bam format protocol for details. - [Total] Fraction of Properly paired // Ratio of properly paired reads against mapped reads - [Total] Read and mate paired // Read (read1) and mate read (read2) paired. - [Total] Fraction of Read and mate paired // Ratio of read and mate paired against mapped reads - [Total] Singletons // Read mapped but mate read unmapped, and vice versa. - [Total] Read and mate map to diff chr // Read and mate read mapped to different chromosome, usually because mapping error and structure variants. - [Total] Read1 // First reads in mate paired sequencing - [Total] Read2 // Mate reads - [Total] Read1(rmdup) // First reads after remove duplications. - [Total] Read2(rmdup) // Mate reads after remove duplications. - [Total] forward strand reads // Number of forward strand reads. - [Total] backward strand reads // Number of backward strand reads. - [Total] PCR duplicate reads // PCR duplications. - [Total] Fraction of PCR duplicate reads // Ratio of PCR duplications. - [Total] Map quality cutoff value // Cutoff map quality score, this value can be set by -q. default is 20, because some variants caller like GATK only consider high quality reads. - [Total] MapQuality above cutoff reads // Number of reads with higher or equal quality score than cutoff value. - [Total] Fraction of MapQ reads in all reads // Ratio of reads with higher or equal Q score against raw reads. - [Total] Fraction of MapQ reads in mapped reads // Ratio of reads with higher or equal Q score against mapped reads. - [Target] Target Reads // Number of reads covered target region (specified by bed file). - [Target] Fraction of Target Reads in all reads // Ratio of target reads against raw reads. - [Target] Fraction of Target Reads in mapped reads // Ratio of target reads against mapped reads. - [Target] Target Data(Mb) // Total bases covered target region. If a read covered target region partly, only the covered bases will be counted. - [Target] Target Data Rmdup(Mb) // Total bases covered target region after remove PCR duplications. - [Target] Fraction of Target Data in all data // Ratio of target bases against raw bases. - [Target] Fraction of Target Data in mapped data // Ratio of target bases against mapped bases. - [Target] Len of region // The length of target regions. - [Target] Average depth // Average depth of target regions. Calculated by "target bases / length of regions". - [Target] Average depth(rmdup) // Average depth of target regions after remove PCR duplications. - [Target] Coverage (>0x) // Ratio of bases with depth greater than 0x in target regions, which also means the ratio of covered regions in target regions. - [Target] Coverage (>=4x) // Ratio of bases with depth greater than or equal to 4x in target regions. - [Target] Coverage (>=10x) // Ratio of bases with depth greater than or equal to 10x in target regions. - [Target] Coverage (>=30x) // Ratio of bases with depth greater than or equal to 30x in target regions. - [Target] Coverage (>=100x) // Ratio of bases with depth greater than or equal to 100x in target regions. - [Target] Coverage (>=Nx) // This is addtional line for user self-defined cutoff value, see --cutoffdepth - [Target] Target Region Count // Number of target regions. In normal practise,it is the total number of exomes. - [Target] Region covered > 0x // The number of these regions with average depth greater than 0x. - [Target] Fraction Region covered > 0x // Ratio of these regions with average depth greater than 0x. - [Target] Fraction Region covered >= 4x // Ratio of these regions with average depth greater than or equal to 4x. - [Target] Fraction Region covered >= 10x // Ratio of these regions with average depth greater than or equal to 10x. - [Target] Fraction Region covered >= 30x // Ratio of these regions with average depth greater than or equal to 30x. - [Target] Fraction Region covered >= 100x // Ratio of these regions with average depth greater than or equal to 100x. - [flank] flank size // The flank size will be count. 200 bp in default. Oligos could also capture the nearby regions of target regions. - [flank] Len of region (not include target region) // The length of flank regions (target regions will not be count). - [flank] Average depth // Average depth of flank regions. - [flank] flank Reads // The total number of reads covered the flank regions. Note: some reads covered the edge of target regions, will be count in flank regions also. - [flank] Fraction of flank Reads in all reads // Ratio of reads covered in flank regions against raw reads. - [flank] Fraction of flank Reads in mapped reads // Ration of reads covered in flank regions against mapped reads. - [flank] flank Data(Mb) // Total bases in the flank regions. - [flank] Fraction of flank Data in all data // Ratio of total bases in the flank regions against raw data. - [flank] Fraction of flank Data in mapped data // Ratio of total bases in the flank regions against mapped data. - [flank] Coverage (>0x) // Ratio of flank bases with depth greater than 0x. - [flank] Coverage (>=4x) // Ratio of flank bases with depth greater than or equal to 4x. - [flank] Coverage (>=10x) // Ratio of flank bases with depth greater than or equal to 10x. - [flank] Coverage (>=30x) // Ratio of flank bases with depth greater than or equal to 30x. - - -**cumu.plot** - -Depth distrbution for plot. - -**insert.plot** - -Inferred insert size distribution for plot. - -**chromosome.report** - -Depth and coverage information of each chromosome. - -**region.tsv.gz** - -For each region in probe file (in.bed), average depth, median - -depth and coverage of these regions will be listed in the file. - -**depth.tsv.gz** - -For each position in the probe file(in.bed), three kinds depth - -value will be calculated, including raw depth, rmdup depth and - -the coverage depth. The raw depth is retrieved from input bam - -file(s) without any restriction. And the rmdup depth is - -calculated after remove duplicated reads and secondary alignment - -reads and low map quality reads(mapQ < 20), this value is similar - -with the output depth of `samtools depth`. The coverage depth is +```bash +# Debian/Ubuntu +apt-get install libhts-dev zlib1g-dev +# macOS (Homebrew) +brew install htslib +``` -the raw depth with consider of deletion region, so this value +## Install -should be equal to or greated than the raw depth. We usw raw depth +```bash +git clone https://github.com/pzweuj/xamdst +cd xamdst +make +``` -to stat the coverage information in the coverage.report file. If +Or use the Docker image: -you want use rmdup depth to calculate the coverage, please use +```bash +docker pull ghcr.io/pzweuj/mapping:2026Jan +``` -the parameter "--use_rmdup". +## Usage -**uncover.bed** +```bash +# basic +xamdst -p probe.bed -o ./ in1.bam -This bed file contains the bad covered or uncovered region in the +# multi-threaded +xamdst -p probe.bed -o ./ --threads 8 in1.bam -input bam file(s) against the probe file. Set the cutoff value of +# CRAM (reference required) +xamdst -p probe.bed -o ./ in1.cram -T ref.fa -uncover by parameter "--uncover" +# streaming from stdin +samtools view in1.bam -u | xamdst -p probe.bed -o ./ - +``` +## Options + +| Option | Default | Description | +| --- | --- | --- | +| `-p, --bed FILE` | — | target / probe regions (merged first). **Required** | +| `-o, --outdir DIR` | — | output directory (created if missing). **Required** | +| `-T, --reference FILE` | — | reference FASTA, **required for CRAM** | +| `-f, --flank N` | 200 | flank bp to stat around each region | +| `-q, --mapthres N` | 20 | mapQ cutoff (≥ counted as quality reads) | +| `--maxdepth N` | 0 | cap depth in the cumulative distribution (0 = no cap) | +| `--cutoffdepth D1,D2,...` | 0 | report coverage at these depths (max 10) | +| `--depthratio R1,R2,...` | 0.2,0.5 | coverage at ratios of average depth | +| `--isize N` | 2000 | insert-size cutoff for the distribution | +| `--uncover N` | 5 | positions below this depth go to `uncover.bed` | +| `--bamout FILE` | — | export target reads to this BAM | +| `--threads N` | 0 | threads for BAM/CRAM I/O (0 = single-threaded) | +| `-1` | off | BED begin position is 1-based | +| `-h, --help` | | show help | +| `-v, --version` | | show version | + +## Output + +Eight files are written to the output directory: + +| File | Content | +| --- | --- | +| `coverage.report` | coverage & read statistics for target and flank regions | +| `coverage.report.json` | same report in JSON | +| `cumu.plot` | depth cumulative distribution | +| `insert.plot` | inferred insert-size distribution | +| `chromosome.report` | per-chromosome coverage | +| `region.tsv.gz` | mean / median depth and coverage per region (bgzip + tabix) | +| `depth.tsv.gz` | per-position raw / rmdup / coverage depth (bgzip + tabix) | +| `uncover.bed` | poorly covered or uncovered regions | + +### `depth.tsv.gz` columns + +| Column | Meaning | +| --- | --- | +| chromosome | chromosome name | +| position | 1-based position | +| raw depth | unfiltered depth | +| rmdup depth | deduplicated, primary, mapQ ≥ cutoff (default 20) — close to `samtools depth` | +| coverage depth | raw depth plus CIGAR deletions (≥ raw depth) | + +Coverage statistics in `coverage.report` use the **coverage depth**. Pass +`--depthratio` / `--cutoffdepth` for custom thresholds. + +## Notes + +- All input files must use the same reference and sort order. +- The first file's header is used for output. +- `region.tsv.gz` and `depth.tsv.gz` can be indexed with `tabix -p bed` / `tabix -p vcf`. + +Homepage: diff --git a/bgzf.c b/bgzf.c deleted file mode 100644 index 9833414..0000000 --- a/bgzf.c +++ /dev/null @@ -1,555 +0,0 @@ -/* The MIT License - - Copyright (c) 2008 Broad Institute / Massachusetts Institute of Technology - 2011 Attractive Chaos - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in - all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN - THE SOFTWARE. -*/ - -#include -#include -#include -#include -#include -#include -#include "bgzf.h" - -#ifdef _USE_KNETFILE -#include "knetfile.h" -typedef knetFile *_bgzf_file_t; -#define _bgzf_open(fn, mode) knet_open(fn, mode) -#define _bgzf_dopen(fp, mode) knet_dopen(fp, mode) -#define _bgzf_close(fp) knet_close(fp) -#define _bgzf_fileno(fp) ((fp)->fd) -#define _bgzf_tell(fp) knet_tell(fp) -#define _bgzf_seek(fp, offset, whence) knet_seek(fp, offset, whence) -#define _bgzf_read(fp, buf, len) knet_read(fp, buf, len) -#define _bgzf_write(fp, buf, len) knet_write(fp, buf, len) -#else // ~defined(_USE_KNETFILE) -#if defined(_WIN32) || defined(_MSC_VER) -#define ftello(fp) ftell(fp) -#define fseeko(fp, offset, whence) fseek(fp, offset, whence) -#else // ~defined(_WIN32) -extern off_t ftello(FILE *stream); -extern int fseeko(FILE *stream, off_t offset, int whence); -#endif // ~defined(_WIN32) -typedef FILE *_bgzf_file_t; -#define _bgzf_open(fn, mode) fopen(fn, mode) -#define _bgzf_dopen(fp, mode) fdopen(fp, mode) -#define _bgzf_close(fp) fclose(fp) -#define _bgzf_fileno(fp) fileno(fp) -#define _bgzf_tell(fp) ftello(fp) -#define _bgzf_seek(fp, offset, whence) fseeko(fp, offset, whence) -#define _bgzf_read(fp, buf, len) fread(buf, 1, len, fp) -#define _bgzf_write(fp, buf, len) fwrite(buf, 1, len, fp) -#endif // ~define(_USE_KNETFILE) - -#define BLOCK_HEADER_LENGTH 18 -#define BLOCK_FOOTER_LENGTH 8 - -/* BGZF/GZIP header (speciallized from RFC 1952; little endian): - +---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+ - | 31|139| 8| 4| 0| 0|255| 6| 66| 67| 2|BLK_LEN| - +---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+ -*/ -static const uint8_t g_magic[19] = "\037\213\010\4\0\0\0\0\0\377\6\0\102\103\2\0\0\0"; - -#ifdef BGZF_CACHE -typedef struct { - int size; - uint8_t *block; - int64_t end_offset; -} cache_t; -#include "khash.h" -KHASH_MAP_INIT_INT64(cache, cache_t) -#endif - -static inline void packInt16(uint8_t *buffer, uint16_t value) -{ - buffer[0] = value; - buffer[1] = value >> 8; -} - -static inline int unpackInt16(const uint8_t *buffer) -{ - return buffer[0] | buffer[1] << 8; -} - -static 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 BGZF *bgzf_read_init() -{ - BGZF *fp; - fp = calloc(1, sizeof(BGZF)); - fp->open_mode = 'r'; - fp->uncompressed_block = malloc(BGZF_MAX_BLOCK_SIZE); - fp->compressed_block = malloc(BGZF_MAX_BLOCK_SIZE); -#ifdef BGZF_CACHE - fp->cache = kh_init(cache); -#endif - return fp; -} - -static BGZF *bgzf_write_init(int compress_level) // compress_level==-1 for the default level -{ - BGZF *fp; - fp = calloc(1, sizeof(BGZF)); - fp->open_mode = 'w'; - fp->uncompressed_block = malloc(BGZF_MAX_BLOCK_SIZE); - fp->compressed_block = malloc(BGZF_MAX_BLOCK_SIZE); - 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; - return fp; -} -// get the compress level from the mode string -static int mode2level(const char *__restrict mode) -{ - 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 compress_level; -} - -BGZF *bgzf_open(const char *path, const char *mode) -{ - BGZF *fp = 0; - if (strchr(mode, 'r') || strchr(mode, 'R')) { - _bgzf_file_t fpr; - if ((fpr = _bgzf_open(path, "r")) == 0) return 0; - fp = bgzf_read_init(); - fp->fp = fpr; - } else if (strchr(mode, 'w') || strchr(mode, 'W')) { - FILE *fpw; - if ((fpw = fopen(path, "w")) == 0) return 0; - fp = bgzf_write_init(mode2level(mode)); - fp->fp = fpw; - } - return fp; -} - -BGZF *bgzf_dopen(int fd, const char *mode) -{ - BGZF *fp = 0; - if (strchr(mode, 'r') || strchr(mode, 'R')) { - _bgzf_file_t fpr; - if ((fpr = _bgzf_dopen(fd, "r")) == 0) return 0; - fp = bgzf_read_init(); - fp->fp = fpr; - } else if (strchr(mode, 'w') || strchr(mode, 'W')) { - FILE *fpw; - if ((fpw = fdopen(fd, "w")) == 0) return 0; - fp = bgzf_write_init(mode2level(mode)); - fp->fp = fpw; - } - return fp; -} - -// Deflate the block in fp->uncompressed_block into fp->compressed_block. Also adds an extra field that stores the compressed block length. -static int deflate_block(BGZF *fp, int block_length) -{ - uint8_t *buffer = fp->compressed_block; - int buffer_size = BGZF_BLOCK_SIZE; - int input_length = block_length; - int compressed_length = 0; - int remaining; - uint32_t crc; - - assert(block_length <= BGZF_BLOCK_SIZE); // guaranteed by the caller - memcpy(buffer, g_magic, BLOCK_HEADER_LENGTH); // the last two bytes are a place holder for the length of the block - while (1) { // loop to retry for blocks that do not compress enough - int status; - 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; - status = deflateInit2(&zs, fp->compress_level, Z_DEFLATED, -15, 8, Z_DEFAULT_STRATEGY); // -15 to disable zlib header/footer - if (status != Z_OK) { - fp->errcode |= BGZF_ERR_ZLIB; - return -1; - } - status = deflate(&zs, Z_FINISH); - if (status != Z_STREAM_END) { // not compressed enough - deflateEnd(&zs); // reset the stream - if (status == Z_OK) { // reduce the size and recompress - input_length -= 1024; - assert(input_length > 0); // logically, this should not happen - continue; - } - fp->errcode |= BGZF_ERR_ZLIB; - return -1; - } - if (deflateEnd(&zs) != Z_OK) { - fp->errcode |= BGZF_ERR_ZLIB; - return -1; - } - compressed_length = zs.total_out; - compressed_length += BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH; - assert(compressed_length <= BGZF_BLOCK_SIZE); - break; - } - - assert(compressed_length > 0); - packInt16((uint8_t*)&buffer[16], compressed_length - 1); // write the compressed_length; -1 to fit 2 bytes - 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); - - remaining = block_length - input_length; - if (remaining > 0) { - assert(remaining <= input_length); - memcpy(fp->uncompressed_block, fp->uncompressed_block + input_length, remaining); - } - fp->block_offset = remaining; - return compressed_length; -} - -// Inflate the block in fp->compressed_block into fp->uncompressed_block -static int inflate_block(BGZF* fp, int block_length) -{ - z_stream zs; - 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 = BGZF_BLOCK_SIZE; - - if (inflateInit2(&zs, -15) != Z_OK) { - fp->errcode |= BGZF_ERR_ZLIB; - return -1; - } - if (inflate(&zs, Z_FINISH) != Z_STREAM_END) { - inflateEnd(&zs); - fp->errcode |= BGZF_ERR_ZLIB; - return -1; - } - if (inflateEnd(&zs) != Z_OK) { - fp->errcode |= BGZF_ERR_ZLIB; - return -1; - } - return zs.total_out; -} - -static int check_header(const uint8_t *header) -{ - return (header[0] == 31 && header[1] == 139 && header[2] == 8 && (header[3] & 4) != 0 - && unpackInt16((uint8_t*)&header[10]) == 6 - && header[12] == 'B' && header[13] == 'C' - && unpackInt16((uint8_t*)&header[14]) == 2); -} - -#ifdef BGZF_CACHE -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, BGZF_BLOCK_SIZE); - _bgzf_seek((_bgzf_file_t)fp->fp, p->end_offset, SEEK_SET); - 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 (BGZF_BLOCK_SIZE >= fp->cache_size) return; - if ((kh_size(h) + 1) * BGZF_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(BGZF_BLOCK_SIZE); - memcpy(kh_val(h, k).block, fp->uncompressed_block, BGZF_BLOCK_SIZE); -} -#else -static void free_cache(BGZF *fp) {} -static int load_block_from_cache(BGZF *fp, int64_t block_address) {return 0;} -static void cache_block(BGZF *fp, int size) {} -#endif - -int bgzf_read_block(BGZF *fp) -{ - uint8_t header[BLOCK_HEADER_LENGTH], *compressed_block; - int count, size = 0, block_length, remaining; - int64_t block_address; - block_address = _bgzf_tell((_bgzf_file_t)fp->fp); - if (load_block_from_cache(fp, block_address)) return 0; - count = _bgzf_read(fp->fp, header, sizeof(header)); - if (count == 0) { // no data read - fp->block_length = 0; - return 0; - } - if (count != sizeof(header) || !check_header(header)) { - fp->errcode |= BGZF_ERR_HEADER; - return -1; - } - size = count; - block_length = unpackInt16((uint8_t*)&header[16]) + 1; // +1 because when writing this number, we used "-1" - compressed_block = (uint8_t*)fp->compressed_block; - memcpy(compressed_block, header, BLOCK_HEADER_LENGTH); - remaining = block_length - BLOCK_HEADER_LENGTH; - count = _bgzf_read(fp->fp, &compressed_block[BLOCK_HEADER_LENGTH], remaining); - if (count != remaining) { - fp->errcode |= BGZF_ERR_IO; - return -1; - } - size += count; - if ((count = inflate_block(fp, block_length)) < 0) return -1; - if (fp->block_length != 0) fp->block_offset = 0; // Do not reset offset if this read follows a seek. - fp->block_address = block_address; - fp->block_length = count; - cache_block(fp, size); - return 0; -} - -ssize_t bgzf_read(BGZF *fp, void *data, ssize_t length) -{ - ssize_t bytes_read = 0; - uint8_t *output = data; - if (length <= 0) return 0; - assert(fp->open_mode == 'r'); - while (bytes_read < length) { - int copy_length, available = fp->block_length - fp->block_offset; - uint8_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 = length - bytes_read < available? 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) { - fp->block_address = _bgzf_tell((_bgzf_file_t)fp->fp); - fp->block_offset = fp->block_length = 0; - } - return bytes_read; -} - -int bgzf_flush(BGZF *fp) -{ - assert(fp->open_mode == 'w'); - while (fp->block_offset > 0) { - int block_length; - block_length = deflate_block(fp, fp->block_offset); - if (block_length < 0) return -1; - if (fwrite(fp->compressed_block, 1, block_length, fp->fp) != block_length) { - fp->errcode |= BGZF_ERR_IO; // possibly truncated file - return -1; - } - fp->block_address += block_length; - } - return 0; -} - -int bgzf_flush_try(BGZF *fp, ssize_t size) -{ - if (fp->block_offset + size > BGZF_BLOCK_SIZE) - return bgzf_flush(fp); - return -1; -} - -ssize_t bgzf_write(BGZF *fp, const void *data, ssize_t length) -{ - const uint8_t *input = data; - int block_length = BGZF_BLOCK_SIZE, bytes_written; - assert(fp->open_mode == 'w'); - input = data; - bytes_written = 0; - while (bytes_written < length) { - uint8_t* buffer = fp->uncompressed_block; - int copy_length = block_length - fp->block_offset < length - bytes_written? block_length - fp->block_offset : length - bytes_written; - 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 && bgzf_flush(fp)) break; - } - return bytes_written; -} - -int bgzf_close(BGZF* fp) -{ - int ret, count, block_length; - if (fp == 0) return -1; - if (fp->open_mode == 'w') { - if (bgzf_flush(fp) != 0) return -1; - block_length = deflate_block(fp, 0); // write an empty block - count = fwrite(fp->compressed_block, 1, block_length, fp->fp); - if (fflush(fp->fp) != 0) { - fp->errcode |= BGZF_ERR_IO; - return -1; - } - } - ret = fp->open_mode == 'w'? fclose(fp->fp) : _bgzf_close(fp->fp); - if (ret != 0) return -1; - 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; - offset = _bgzf_tell((_bgzf_file_t)fp->fp); - if (_bgzf_seek(fp->fp, -28, SEEK_END) < 0) return 0; - _bgzf_read(fp->fp, buf, 28); - _bgzf_seek(fp->fp, offset, SEEK_SET); - 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' || where != SEEK_SET) { - fp->errcode |= BGZF_ERR_MISUSE; - return -1; - } - block_offset = pos & 0xFFFF; - block_address = pos >> 16; - if (_bgzf_seek(fp->fp, block_address, SEEK_SET) < 0) { - fp->errcode |= BGZF_ERR_IO; - return -1; - } - fp->block_length = 0; // indicates current block has not been loaded - fp->block_address = block_address; - fp->block_offset = block_offset; - return 0; -} - -int bgzf_is_bgzf(const char *fn) -{ - uint8_t buf[16]; - int n; - _bgzf_file_t fp; - if ((fp = _bgzf_open(fn, "r")) == 0) return 0; - n = _bgzf_read(fp, buf, 16); - _bgzf_close(fp); - if (n != 16) return 0; - return memcmp(g_magic, buf, 16) == 0? 1 : 0; -} - -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) { - fp->block_address = _bgzf_tell((_bgzf_file_t)fp->fp); - fp->block_offset = 0; - fp->block_length = 0; - } - return c; -} - -#ifndef kroundup32 -#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) -#endif - -int bgzf_getline(BGZF *fp, int delim, 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] != delim; ++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) { - fp->block_address = _bgzf_tell((_bgzf_file_t)fp->fp); - 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; -} diff --git a/bgzf.h b/bgzf.h deleted file mode 100644 index 29fe0e5..0000000 --- a/bgzf.h +++ /dev/null @@ -1,196 +0,0 @@ -/* The MIT License - - Copyright (c) 2008 Broad Institute / Massachusetts Institute of Technology - 2011 Attractive Chaos - - 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. -*/ - -/* The BGZF library was originally written by Bob Handsaker from the Broad - * Institute. It was later improved by the SAMtools developers. */ - -#ifndef __BGZF_H -#define __BGZF_H - -#include -#include -#include - -#define BGZF_BLOCK_SIZE 0x10000 -#define BGZF_MAX_BLOCK_SIZE 0x10000 - -#define BGZF_ERR_ZLIB 1 -#define BGZF_ERR_HEADER 2 -#define BGZF_ERR_IO 4 -#define BGZF_ERR_MISUSE 8 - -typedef struct { - int open_mode:8, compress_level:8, errcode:16; - int cache_size; - int block_length, block_offset; - int64_t block_address; - void *uncompressed_block, *compressed_block; - void *cache; // a pointer to a hash table - void *fp; // actual file handler; FILE* on writing; FILE* or knetFile* on reading -} BGZF; - -#ifndef KSTRING_T -#define KSTRING_T kstring_t -typedef struct __kstring_t { - size_t l, m; - char *s; -} kstring_t; -#endif - -#ifdef __cplusplus -extern "C" { -#endif - - /****************** - * Basic routines * - ******************/ - - /** - * Open an existing file descriptor for reading or writing. - * - * @param fd file descriptor - * @param mode mode matching /[rwu0-9]+/: 'r' for reading, 'w' for writing and a digit specifies - * the zlib compression level; if both 'r' and 'w' are present, 'w' is ignored. - * @return BGZF file handler; 0 on error - */ - BGZF* bgzf_dopen(int fd, const char *mode); - - #define bgzf_fdopen(fd, mode) bgzf_dopen((fd), (mode)) // for backward compatibility - - /** - * Open the specified file for reading or writing. - */ - BGZF* bgzf_open(const char* path, const char *mode); - - /** - * Close the BGZF and free all associated resources. - * - * @param fp BGZF file handler - * @return 0 on success and -1 on error - */ - int bgzf_close(BGZF *fp); - - /** - * Read up to _length_ bytes from the file storing into _data_. - * - * @param fp BGZF file handler - * @param data data array to read into - * @param length size of data to read - * @return number of bytes actually read; 0 on end-of-file and -1 on error - */ - ssize_t bgzf_read(BGZF *fp, void *data, ssize_t length); - - /** - * Write _length_ bytes from _data_ to the file. - * - * @param fp BGZF file handler - * @param data data array to write - * @param length size of data to write - * @return number of bytes actually written; -1 on error - */ - ssize_t bgzf_write(BGZF *fp, const void *data, ssize_t length); - - /** - * Write the data in the buffer to the file. - */ - int bgzf_flush(BGZF *fp); - - /** - * 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. - */ - #define bgzf_tell(fp) ((fp->block_address << 16) | (fp->block_offset & 0xFFFF)) - - /** - * Set the file to read from the location specified by _pos_. - * - * @param fp BGZF file handler - * @param pos virtual file offset returned by bgzf_tell() - * @param whence must be SEEK_SET - * @return 0 on success and -1 on error - */ - int64_t bgzf_seek(BGZF *fp, int64_t pos, int whence); - - /** - * Check if the BGZF end-of-file (EOF) marker is present - * - * @param fp BGZF file handler opened for reading - * @return 1 if EOF is present; 0 if not or on I/O error - */ - int bgzf_check_EOF(BGZF *fp); - - /** - * Check if a file is in the BGZF format - * - * @param fn file name - * @return 1 if _fn_ is BGZF; 0 if not or on I/O error - */ - int bgzf_is_bgzf(const char *fn); - - /********************* - * Advanced routines * - *********************/ - - /** - * Set the cache size. Only effective when compiled with -DBGZF_CACHE. - * - * @param fp BGZF file handler - * @param size size of cache in bytes; 0 to disable caching (default) - */ - void bgzf_set_cache_size(BGZF *fp, int size); - - /** - * Flush the file if the remaining buffer size is smaller than _size_ - */ - int bgzf_flush_try(BGZF *fp, ssize_t size); - - /** - * Read one byte from a BGZF file. It is faster than bgzf_read() - * @param fp BGZF file handler - * @return byte read; -1 on end-of-file or error - */ - int bgzf_getc(BGZF *fp); - - /** - * Read one line from a BGZF file. It is faster than bgzf_getc() - * - * @param fp BGZF file handler - * @param delim delimitor - * @param str string to write to; must be initialized - * @return length of the string; 0 on end-of-file; negative on error - */ - int bgzf_getline(BGZF *fp, int delim, kstring_t *str); - - /** - * Read the next BGZF block. - */ - int bgzf_read_block(BGZF *fp); - -#ifdef __cplusplus -} -#endif - -#endif diff --git a/count.h b/count.h index fb21eac..aa18473 100644 --- a/count.h +++ b/count.h @@ -51,17 +51,21 @@ count64_t; #define count32_init(c) do { \ (c) = (count32_t *)malloc(sizeof(count32_t)); \ + if ((c) == NULL) errabort("[count32_init] out of memory"); \ (c)->m = 0; \ (c)->n = 64; \ (c)->a = malloc(sizeof(uint32_t)*(c)->n); \ + if ((c)->a == NULL) errabort("[count32_init] out of memory"); \ memset((c)->a, 0, (c)->n *sizeof(uint32_t)); \ } while(0) #define count64_init(c) do { \ (c) = (count64_t*)malloc(sizeof(count64_t)); \ + if ((c) == NULL) errabort("[count64_init] out of memory"); \ (c)->m = 0; \ (c)->n = 64; \ (c)->a = malloc(sizeof(uint64_t)*(c)->n); \ + if ((c)->a == NULL) errabort("[count64_init] out of memory"); \ memset((c)->a, 0, (c)->n *sizeof(uint64_t)); \ } while(0) @@ -78,13 +82,17 @@ count64_t; if ((c)->n <= (c)->m) \ { \ (c)->n = (c)->m + 1024; \ - (c)->a = realloc((c)->a, sizeof(type) * (c)->n); \ + type *_r = (type*)realloc((c)->a, sizeof(type) * (c)->n); \ + if (_r == NULL) errabort("[count_increase] out of memory"); \ + (c)->a = _r; \ memset((c)->a + (c)->m, 0, ((c)->n-(c)->m)*sizeof(type)); \ } \ if((c)->n <= d+1) \ { \ (c)->n = d+1024; \ - (c)->a = realloc((c)->a, sizeof(type) *(c)->n); \ + type *_r = (type*)realloc((c)->a, sizeof(type) *(c)->n); \ + if (_r == NULL) errabort("[count_increase] out of memory"); \ + (c)->a = _r; \ memset((c)->a+(c)->m, 0, ((c)->n - (c)->m)*sizeof(type)); \ (c)->m = d+1; \ } \ @@ -95,18 +103,22 @@ count64_t; #define count_resize(c, l, type) do { \ if((c)->n < l) \ { \ - (c)->a = realloc((c)->a, sizeof(type) * l); \ + type *_r = (type*)realloc((c)->a, sizeof(type) * l); \ + if (_r == NULL) errabort("[count_resize] out of memory"); \ + (c)->a = _r; \ int i; \ for (i = (c)->n; i < l; ++i) (c)->a[i] = 0; \ (c)->n = l; \ } \ - } while(0); + } while(0) #define count_increaseN(c, d, cnt, type) do { \ if((c)->n <= d+1) \ { \ (c)->n = d+1024; \ - (c)->a = realloc((c)->a, sizeof(type) *(c)->n); \ + type *_r = (type*)realloc((c)->a, sizeof(type) *(c)->n); \ + if (_r == NULL) errabort("[count_increaseN] out of memory"); \ + (c)->a = _r; \ memset((c)->a+(c)->m, 0, ((c)->n - (c)->m)* sizeof(type)); \ (c)->m = d+1; \ } \ diff --git a/samlib/bam.c b/samlib/bam.c deleted file mode 100644 index b00d6a6..0000000 --- a/samlib/bam.c +++ /dev/null @@ -1,474 +0,0 @@ -#include -#include -#include -#include -#include "bam.h" -#include "bam_endian.h" -#include "kstring.h" -#include "sam_header.h" - -int bam_is_be = 0, bam_verbose = 2, bam_no_B = 0; -char *bam_flag2char_table = "pPuUrR12sfd\0\0\0\0\0"; - -/************************** - * CIGAR related routines * - **************************/ - -uint32_t bam_calend(const bam1_core_t *c, const uint32_t *cigar) -{ - int k, end = c->pos; - for (k = 0; k < c->n_cigar; ++k) { - int op = bam_cigar_op(cigar[k]); - int len = bam_cigar_oplen(cigar[k]); - if (op == BAM_CBACK) { // move backward - int l, u, v; - if (k == c->n_cigar - 1) break; // skip trailing 'B' - for (l = k - 1, u = v = 0; l >= 0; --l) { - int op1 = bam_cigar_op(cigar[l]); - int len1 = bam_cigar_oplen(cigar[l]); - if (bam_cigar_type(op1)&1) { // consume query - if (u + len1 >= len) { // stop - if (bam_cigar_type(op1)&2) v += len - u; - break; - } else u += len1; - } - if (bam_cigar_type(op1)&2) v += len1; - } - end = l < 0? c->pos : end - v; - } else if (bam_cigar_type(op)&2) end += bam_cigar_oplen(cigar[k]); - } - return end; -} - -int32_t bam_cigar2qlen(const bam1_core_t *c, const uint32_t *cigar) -{ - uint32_t k; - int32_t l = 0; - for (k = 0; k < c->n_cigar; ++k) - if (bam_cigar_type(bam_cigar_op(cigar[k]))&1) - l += bam_cigar_oplen(cigar[k]); - return l; -} - -/******************** - * BAM I/O routines * - ********************/ - -bam_header_t *bam_header_init() -{ - bam_is_be = bam_is_big_endian(); - return (bam_header_t*)calloc(1, sizeof(bam_header_t)); -} - -void bam_header_destroy(bam_header_t *header) -{ - int32_t i; - extern void bam_destroy_header_hash(bam_header_t *header); - if (header == 0) return; - if (header->target_name) { - for (i = 0; i < header->n_targets; ++i) - free(header->target_name[i]); - free(header->target_name); - free(header->target_len); - } - free(header->text); - if (header->dict) sam_header_free(header->dict); - if (header->rg2lib) sam_tbl_destroy(header->rg2lib); - bam_destroy_header_hash(header); - free(header); -} - -bam_header_t *bam_header_read(bamFile fp) -{ - bam_header_t *header; - char buf[4]; - int magic_len; - int32_t i = 1, name_len; - // check EOF - i = bgzf_check_EOF(fp); - if (i < 0) { - // If the file is a pipe, checking the EOF marker will *always* fail - // with ESPIPE. Suppress the error message in this case. - if (errno != ESPIPE) perror("[bam_header_read] bgzf_check_EOF"); - } - else if (i == 0) fprintf(stderr, "[bam_header_read] EOF marker is absent. The input is probably truncated.\n"); - // read "BAM1" - magic_len = bam_read(fp, buf, 4); - if (magic_len != 4 || strncmp(buf, "BAM\001", 4) != 0) { - fprintf(stderr, "[bam_header_read] invalid BAM binary header (this is not a BAM file).\n"); - return 0; - } - header = bam_header_init(); - // read plain text and the number of reference sequences - bam_read(fp, &header->l_text, 4); - if (bam_is_be) bam_swap_endian_4p(&header->l_text); - header->text = (char*)calloc(header->l_text + 1, 1); - bam_read(fp, header->text, header->l_text); - bam_read(fp, &header->n_targets, 4); - if (bam_is_be) bam_swap_endian_4p(&header->n_targets); - // read reference sequence names and lengths - header->target_name = (char**)calloc(header->n_targets, sizeof(char*)); - header->target_len = (uint32_t*)calloc(header->n_targets, 4); - for (i = 0; i != header->n_targets; ++i) { - bam_read(fp, &name_len, 4); - if (bam_is_be) bam_swap_endian_4p(&name_len); - header->target_name[i] = (char*)calloc(name_len, 1); - bam_read(fp, header->target_name[i], name_len); - bam_read(fp, &header->target_len[i], 4); - if (bam_is_be) bam_swap_endian_4p(&header->target_len[i]); - } - return header; -} - -int bam_header_write(bamFile fp, const bam_header_t *header) -{ - char buf[4]; - int32_t i, name_len, x; - // write "BAM1" - strncpy(buf, "BAM\001", 4); - bam_write(fp, buf, 4); - // write plain text and the number of reference sequences - if (bam_is_be) { - x = bam_swap_endian_4(header->l_text); - bam_write(fp, &x, 4); - if (header->l_text) bam_write(fp, header->text, header->l_text); - x = bam_swap_endian_4(header->n_targets); - bam_write(fp, &x, 4); - } else { - bam_write(fp, &header->l_text, 4); - if (header->l_text) bam_write(fp, header->text, header->l_text); - bam_write(fp, &header->n_targets, 4); - } - // write sequence names and lengths - for (i = 0; i != header->n_targets; ++i) { - char *p = header->target_name[i]; - name_len = strlen(p) + 1; - if (bam_is_be) { - x = bam_swap_endian_4(name_len); - bam_write(fp, &x, 4); - } else bam_write(fp, &name_len, 4); - bam_write(fp, p, name_len); - if (bam_is_be) { - x = bam_swap_endian_4(header->target_len[i]); - bam_write(fp, &x, 4); - } else bam_write(fp, &header->target_len[i], 4); - } - bgzf_flush(fp); - return 0; -} - -static void swap_endian_data(const bam1_core_t *c, int data_len, uint8_t *data) -{ - uint8_t *s; - uint32_t i, *cigar = (uint32_t*)(data + c->l_qname); - s = data + c->n_cigar*4 + c->l_qname + c->l_qseq + (c->l_qseq + 1)/2; - for (i = 0; i < c->n_cigar; ++i) bam_swap_endian_4p(&cigar[i]); - while (s < data + data_len) { - uint8_t type; - s += 2; // skip key - type = toupper(*s); ++s; // skip type - if (type == 'C' || type == 'A') ++s; - else if (type == 'S') { bam_swap_endian_2p(s); s += 2; } - else if (type == 'I' || type == 'F') { bam_swap_endian_4p(s); s += 4; } - else if (type == 'D') { bam_swap_endian_8p(s); s += 8; } - else if (type == 'Z' || type == 'H') { while (*s) ++s; ++s; } - else if (type == 'B') { - int32_t n, Bsize = bam_aux_type2size(*s); - memcpy(&n, s + 1, 4); - if (1 == Bsize) { - } else if (2 == Bsize) { - for (i = 0; i < n; i += 2) - bam_swap_endian_2p(s + 5 + i); - } else if (4 == Bsize) { - for (i = 0; i < n; i += 4) - bam_swap_endian_4p(s + 5 + i); - } - bam_swap_endian_4p(s+1); - } - } -} - -int bam_read1(bamFile fp, bam1_t *b) -{ - bam1_core_t *c = &b->core; - int32_t block_len, ret, i; - uint32_t x[8]; - - assert(BAM_CORE_SIZE == 32); - if ((ret = bam_read(fp, &block_len, 4)) != 4) { - if (ret == 0) return -1; // normal end-of-file - else return -2; // truncated - } - if (bam_read(fp, x, BAM_CORE_SIZE) != BAM_CORE_SIZE) return -3; - if (bam_is_be) { - bam_swap_endian_4p(&block_len); - for (i = 0; i < 8; ++i) bam_swap_endian_4p(x + i); - } - c->tid = x[0]; c->pos = x[1]; - c->bin = x[2]>>16; c->qual = x[2]>>8&0xff; c->l_qname = x[2]&0xff; - c->flag = x[3]>>16; c->n_cigar = x[3]&0xffff; - c->l_qseq = x[4]; - c->mtid = x[5]; c->mpos = x[6]; c->isize = x[7]; - b->data_len = block_len - BAM_CORE_SIZE; - if (b->m_data < b->data_len) { - b->m_data = b->data_len; - kroundup32(b->m_data); - b->data = (uint8_t*)realloc(b->data, b->m_data); - } - if (bam_read(fp, b->data, b->data_len) != b->data_len) return -4; - b->l_aux = b->data_len - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2; - if (bam_is_be) swap_endian_data(c, b->data_len, b->data); - if (bam_no_B) bam_remove_B(b); - return 4 + block_len; -} - -inline int bam_write1_core(bamFile fp, const bam1_core_t *c, int data_len, uint8_t *data) -{ - uint32_t x[8], block_len = data_len + BAM_CORE_SIZE, y; - int i; - assert(BAM_CORE_SIZE == 32); - x[0] = c->tid; - x[1] = c->pos; - x[2] = (uint32_t)c->bin<<16 | c->qual<<8 | c->l_qname; - x[3] = (uint32_t)c->flag<<16 | c->n_cigar; - x[4] = c->l_qseq; - x[5] = c->mtid; - x[6] = c->mpos; - x[7] = c->isize; - bgzf_flush_try(fp, 4 + block_len); - if (bam_is_be) { - for (i = 0; i < 8; ++i) bam_swap_endian_4p(x + i); - y = block_len; - bam_write(fp, bam_swap_endian_4p(&y), 4); - swap_endian_data(c, data_len, data); - } else bam_write(fp, &block_len, 4); - bam_write(fp, x, BAM_CORE_SIZE); - bam_write(fp, data, data_len); - if (bam_is_be) swap_endian_data(c, data_len, data); - return 4 + block_len; -} - -int bam_write1(bamFile fp, const bam1_t *b) -{ - return bam_write1_core(fp, &b->core, b->data_len, b->data); -} - -char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int of) -{ - uint8_t *s = bam1_seq(b), *t = bam1_qual(b); - int i; - const bam1_core_t *c = &b->core; - kstring_t str; - str.l = str.m = 0; str.s = 0; - - kputsn(bam1_qname(b), c->l_qname-1, &str); kputc('\t', &str); - if (of == BAM_OFDEC) { kputw(c->flag, &str); kputc('\t', &str); } - else if (of == BAM_OFHEX) ksprintf(&str, "0x%x\t", c->flag); - else { // BAM_OFSTR - for (i = 0; i < 16; ++i) - if ((c->flag & 1<tid < 0) kputsn("*\t", 2, &str); - else { - if (header) kputs(header->target_name[c->tid] , &str); - else kputw(c->tid, &str); - kputc('\t', &str); - } - kputw(c->pos + 1, &str); kputc('\t', &str); kputw(c->qual, &str); kputc('\t', &str); - if (c->n_cigar == 0) kputc('*', &str); - else { - uint32_t *cigar = bam1_cigar(b); - for (i = 0; i < c->n_cigar; ++i) { - kputw(bam1_cigar(b)[i]>>BAM_CIGAR_SHIFT, &str); - kputc(bam_cigar_opchr(cigar[i]), &str); - } - } - kputc('\t', &str); - if (c->mtid < 0) kputsn("*\t", 2, &str); - else if (c->mtid == c->tid) kputsn("=\t", 2, &str); - else { - if (header) kputs(header->target_name[c->mtid], &str); - else kputw(c->mtid, &str); - kputc('\t', &str); - } - kputw(c->mpos + 1, &str); kputc('\t', &str); kputw(c->isize, &str); kputc('\t', &str); - if (c->l_qseq) { - for (i = 0; i < c->l_qseq; ++i) kputc(bam_nt16_rev_table[bam1_seqi(s, i)], &str); - kputc('\t', &str); - if (t[0] == 0xff) kputc('*', &str); - else for (i = 0; i < c->l_qseq; ++i) kputc(t[i] + 33, &str); - } else kputsn("*\t*", 3, &str); - s = bam1_aux(b); - while (s < b->data + b->data_len) { - uint8_t type, key[2]; - key[0] = s[0]; key[1] = s[1]; - s += 2; type = *s; ++s; - kputc('\t', &str); kputsn((char*)key, 2, &str); kputc(':', &str); - if (type == 'A') { kputsn("A:", 2, &str); kputc(*s, &str); ++s; } - else if (type == 'C') { kputsn("i:", 2, &str); kputw(*s, &str); ++s; } - else if (type == 'c') { kputsn("i:", 2, &str); kputw(*(int8_t*)s, &str); ++s; } - else if (type == 'S') { kputsn("i:", 2, &str); kputw(*(uint16_t*)s, &str); s += 2; } - else if (type == 's') { kputsn("i:", 2, &str); kputw(*(int16_t*)s, &str); s += 2; } - else if (type == 'I') { kputsn("i:", 2, &str); kputuw(*(uint32_t*)s, &str); s += 4; } - else if (type == 'i') { kputsn("i:", 2, &str); kputw(*(int32_t*)s, &str); s += 4; } - else if (type == 'f') { ksprintf(&str, "f:%g", *(float*)s); s += 4; } - else if (type == 'd') { ksprintf(&str, "d:%lg", *(double*)s); s += 8; } - else if (type == 'Z' || type == 'H') { kputc(type, &str); kputc(':', &str); while (*s) kputc(*s++, &str); ++s; } - else if (type == 'B') { - uint8_t sub_type = *(s++); - int32_t n; - memcpy(&n, s, 4); - s += 4; // no point to the start of the array - kputc(type, &str); kputc(':', &str); kputc(sub_type, &str); // write the typing - for (i = 0; i < n; ++i) { - kputc(',', &str); - if ('c' == sub_type || 'c' == sub_type) { kputw(*(int8_t*)s, &str); ++s; } - else if ('C' == sub_type) { kputw(*(uint8_t*)s, &str); ++s; } - else if ('s' == sub_type) { kputw(*(int16_t*)s, &str); s += 2; } - else if ('S' == sub_type) { kputw(*(uint16_t*)s, &str); s += 2; } - else if ('i' == sub_type) { kputw(*(int32_t*)s, &str); s += 4; } - else if ('I' == sub_type) { kputuw(*(uint32_t*)s, &str); s += 4; } - else if ('f' == sub_type) { ksprintf(&str, "%g", *(float*)s); s += 4; } - } - } - } - return str.s; -} - -char *bam_format1(const bam_header_t *header, const bam1_t *b) -{ - return bam_format1_core(header, b, BAM_OFDEC); -} - -void bam_view1(const bam_header_t *header, const bam1_t *b) -{ - char *s = bam_format1(header, b); - puts(s); - free(s); -} - -int bam_validate1(const bam_header_t *header, const bam1_t *b) -{ - char *s; - - if (b->core.tid < -1 || b->core.mtid < -1) return 0; - if (header && (b->core.tid >= header->n_targets || b->core.mtid >= header->n_targets)) return 0; - - if (b->data_len < b->core.l_qname) return 0; - s = memchr(bam1_qname(b), '\0', b->core.l_qname); - if (s != &bam1_qname(b)[b->core.l_qname-1]) return 0; - - // FIXME: Other fields could also be checked, especially the auxiliary data - - return 1; -} - -// FIXME: we should also check the LB tag associated with each alignment -const char *bam_get_library(bam_header_t *h, const bam1_t *b) -{ - const uint8_t *rg; - if (h->dict == 0) h->dict = sam_header_parse2(h->text); - if (h->rg2lib == 0) h->rg2lib = sam_header2tbl(h->dict, "RG", "ID", "LB"); - rg = bam_aux_get(b, "RG"); - return (rg == 0)? 0 : sam_tbl_get(h->rg2lib, (const char*)(rg + 1)); -} - -/************ - * Remove B * - ************/ - -int bam_remove_B(bam1_t *b) -{ - int i, j, end_j, k, l, no_qual; - uint32_t *cigar, *new_cigar; - uint8_t *seq, *qual, *p; - // test if removal is necessary - if (b->core.flag & BAM_FUNMAP) return 0; // unmapped; do nothing - cigar = bam1_cigar(b); - for (k = 0; k < b->core.n_cigar; ++k) - if (bam_cigar_op(cigar[k]) == BAM_CBACK) break; - if (k == b->core.n_cigar) return 0; // no 'B' - if (bam_cigar_op(cigar[0]) == BAM_CBACK) goto rmB_err; // cannot be removed - // allocate memory for the new CIGAR - if (b->data_len + (b->core.n_cigar + 1) * 4 > b->m_data) { // not enough memory - b->m_data = b->data_len + b->core.n_cigar * 4; - kroundup32(b->m_data); - b->data = (uint8_t*)realloc(b->data, b->m_data); - cigar = bam1_cigar(b); // after realloc, cigar may be changed - } - new_cigar = (uint32_t*)(b->data + (b->m_data - b->core.n_cigar * 4)); // from the end of b->data - // the core loop - seq = bam1_seq(b); qual = bam1_qual(b); - no_qual = (qual[0] == 0xff); // test whether base quality is available - i = j = 0; end_j = -1; - for (k = l = 0; k < b->core.n_cigar; ++k) { - int op = bam_cigar_op(cigar[k]); - int len = bam_cigar_oplen(cigar[k]); - if (op == BAM_CBACK) { // the backward operation - int t, u; - if (k == b->core.n_cigar - 1) break; // ignore 'B' at the end of CIGAR - if (len > j) goto rmB_err; // an excessively long backward - for (t = l - 1, u = 0; t >= 0; --t) { // look back - int op1 = bam_cigar_op(new_cigar[t]); - int len1 = bam_cigar_oplen(new_cigar[t]); - if (bam_cigar_type(op1)&1) { // consume the query - if (u + len1 >= len) { // stop - new_cigar[t] -= (len - u) << BAM_CIGAR_SHIFT; - break; - } else u += len1; - } - } - if (bam_cigar_oplen(new_cigar[t]) == 0) --t; // squeeze out the zero-length operation - l = t + 1; - end_j = j; j -= len; - } else { // other CIGAR operations - new_cigar[l++] = cigar[k]; - if (bam_cigar_type(op)&1) { // consume the query - if (i != j) { // no need to copy if i == j - int u, c, c0; - for (u = 0; u < len; ++u) { // construct the consensus - c = bam1_seqi(seq, i+u); - if (j + u < end_j) { // in an overlap - c0 = bam1_seqi(seq, j+u); - if (c != c0) { // a mismatch; choose the better base - if (qual[j+u] < qual[i+u]) { // the base in the 2nd segment is better - bam1_seq_seti(seq, j+u, c); - qual[j+u] = qual[i+u] - qual[j+u]; - } else qual[j+u] -= qual[i+u]; // the 1st is better; reduce base quality - } else qual[j+u] = qual[j+u] > qual[i+u]? qual[j+u] : qual[i+u]; - } else { // not in an overlap; copy over - bam1_seq_seti(seq, j+u, c); - qual[j+u] = qual[i+u]; - } - } - } - i += len, j += len; - } - } - } - if (no_qual) qual[0] = 0xff; // in very rare cases, this may be modified - // merge adjacent operations if possible - for (k = 1; k < l; ++k) - if (bam_cigar_op(new_cigar[k]) == bam_cigar_op(new_cigar[k-1])) - new_cigar[k] += new_cigar[k-1] >> BAM_CIGAR_SHIFT << BAM_CIGAR_SHIFT, new_cigar[k-1] &= 0xf; - // kill zero length operations - for (k = i = 0; k < l; ++k) - if (new_cigar[k] >> BAM_CIGAR_SHIFT) - new_cigar[i++] = new_cigar[k]; - l = i; - // update b - memcpy(cigar, new_cigar, l * 4); // set CIGAR - p = b->data + b->core.l_qname + l * 4; - memmove(p, seq, (j+1)>>1); p += (j+1)>>1; // set SEQ - memmove(p, qual, j); p += j; // set QUAL - memmove(p, bam1_aux(b), b->l_aux); p += b->l_aux; // set optional fields - b->core.n_cigar = l, b->core.l_qseq = j; // update CIGAR length and query length - b->data_len = p - b->data; // update record length - return 0; - -rmB_err: - b->core.flag |= BAM_FUNMAP; - return -1; -} diff --git a/samlib/bam.h b/samlib/bam.h deleted file mode 100644 index b9a1606..0000000 --- a/samlib/bam.h +++ /dev/null @@ -1,795 +0,0 @@ -/* The MIT License - - Copyright (c) 2008-2010 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 */ - -#ifndef BAM_BAM_H -#define BAM_BAM_H - -/*! - @header - - BAM library provides I/O and various operations on manipulating files - in the BAM (Binary Alignment/Mapping) or SAM (Sequence Alignment/Map) - format. It now supports importing from or exporting to SAM, sorting, - merging, generating pileup, and quickly retrieval of reads overlapped - with a specified region. - - @copyright Genome Research Ltd. - */ - -#define BAM_VERSION "0.1.19-44428cd" - -#include -#include -#include -#include - -#ifndef BAM_LITE -#define BAM_VIRTUAL_OFFSET16 -#include "bgzf.h" -/*! @abstract BAM file handler */ -typedef BGZF *bamFile; -#define bam_open(fn, mode) bgzf_open(fn, mode) -#define bam_dopen(fd, mode) bgzf_fdopen(fd, mode) -#define bam_close(fp) bgzf_close(fp) -#define bam_read(fp, buf, size) bgzf_read(fp, buf, size) -#define bam_write(fp, buf, size) bgzf_write(fp, buf, size) -#define bam_tell(fp) bgzf_tell(fp) -#define bam_seek(fp, pos, dir) bgzf_seek(fp, pos, dir) -#else -#define BAM_TRUE_OFFSET -#include -typedef gzFile bamFile; -#define bam_open(fn, mode) gzopen(fn, mode) -#define bam_dopen(fd, mode) gzdopen(fd, mode) -#define bam_close(fp) gzclose(fp) -#define bam_read(fp, buf, size) gzread(fp, buf, size) -/* no bam_write/bam_tell/bam_seek() here */ -#endif - -/*! @typedef - @abstract Structure for the alignment header. - @field n_targets number of reference sequences - @field target_name names of the reference sequences - @field target_len lengths of the referene sequences - @field dict header dictionary - @field hash hash table for fast name lookup - @field rg2lib hash table for @RG-ID -> LB lookup - @field l_text length of the plain text in the header - @field text plain text - - @discussion Field hash points to null by default. It is a private - member. - */ -typedef struct { - int32_t n_targets; - char **target_name; - uint32_t *target_len; - void *dict, *hash, *rg2lib; - uint32_t l_text, n_text; - char *text; -} bam_header_t; - -/*! @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */ -#define BAM_FPAIRED 1 -/*! @abstract the read is mapped in a proper pair */ -#define BAM_FPROPER_PAIR 2 -/*! @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */ -#define BAM_FUNMAP 4 -/*! @abstract the mate is unmapped */ -#define BAM_FMUNMAP 8 -/*! @abstract the read is mapped to the reverse strand */ -#define BAM_FREVERSE 16 -/*! @abstract the mate is mapped to the reverse strand */ -#define BAM_FMREVERSE 32 -/*! @abstract this is read1 */ -#define BAM_FREAD1 64 -/*! @abstract this is read2 */ -#define BAM_FREAD2 128 -/*! @abstract not primary alignment */ -#define BAM_FSECONDARY 256 -/*! @abstract QC failure */ -#define BAM_FQCFAIL 512 -/*! @abstract optical or PCR duplicate */ -#define BAM_FDUP 1024 -// updated 20181031 -#define BAM_FSUPPLEMENTARY 2048 - -#define BAM_OFDEC 0 -#define BAM_OFHEX 1 -#define BAM_OFSTR 2 - -/*! @abstract defautl mask for pileup */ -#define BAM_DEF_MASK (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP) - -#define BAM_CORE_SIZE sizeof(bam1_core_t) - -/** - * Describing how CIGAR operation/length is packed in a 32-bit integer. - */ -#define BAM_CIGAR_SHIFT 4 -#define BAM_CIGAR_MASK ((1 << BAM_CIGAR_SHIFT) - 1) - -/* - CIGAR operations. - */ -/*! @abstract CIGAR: M = match or mismatch*/ -#define BAM_CMATCH 0 -/*! @abstract CIGAR: I = insertion to the reference */ -#define BAM_CINS 1 -/*! @abstract CIGAR: D = deletion from the reference */ -#define BAM_CDEL 2 -/*! @abstract CIGAR: N = skip on the reference (e.g. spliced alignment) */ -#define BAM_CREF_SKIP 3 -/*! @abstract CIGAR: S = clip on the read with clipped sequence - present in qseq */ -#define BAM_CSOFT_CLIP 4 -/*! @abstract CIGAR: H = clip on the read with clipped sequence trimmed off */ -#define BAM_CHARD_CLIP 5 -/*! @abstract CIGAR: P = padding */ -#define BAM_CPAD 6 -/*! @abstract CIGAR: equals = match */ -#define BAM_CEQUAL 7 -/*! @abstract CIGAR: X = mismatch */ -#define BAM_CDIFF 8 -#define BAM_CBACK 9 - -#define BAM_CIGAR_STR "MIDNSHP=XB" -#define BAM_CIGAR_TYPE 0x3C1A7 - -#define bam_cigar_op(c) ((c)&BAM_CIGAR_MASK) -#define bam_cigar_oplen(c) ((c)>>BAM_CIGAR_SHIFT) -#define bam_cigar_opchr(c) (BAM_CIGAR_STR[bam_cigar_op(c)]) -#define bam_cigar_gen(l, o) ((l)<>((o)<<1)&3) // bit 1: consume query; bit 2: consume reference - -/*! @typedef - @abstract Structure for core alignment information. - @field tid chromosome ID, defined by bam_header_t - @field pos 0-based leftmost coordinate - @field bin bin calculated by bam_reg2bin() - @field qual mapping quality - @field l_qname length of the query name - @field flag bitwise flag - @field n_cigar number of CIGAR operations - @field l_qseq length of the query sequence (read) - */ -typedef struct { - int32_t tid; - int32_t pos; - uint32_t bin:16, qual:8, l_qname:8; - uint32_t flag:16, n_cigar:16; - int32_t l_qseq; - int32_t mtid; - int32_t mpos; - int32_t isize; -} bam1_core_t; - -/*! @typedef - @abstract Structure for one alignment. - @field core core information about the alignment - @field l_aux length of auxiliary data - @field data_len current length of bam1_t::data - @field m_data maximum length of bam1_t::data - @field data all variable-length data, concatenated; structure: qname-cigar-seq-qual-aux - - @discussion Notes: - - 1. qname is zero tailing and core.l_qname includes the tailing '\0'. - 2. l_qseq is calculated from the total length of an alignment block - on reading or from CIGAR. - 3. cigar data is encoded 4 bytes per CIGAR operation. - 4. seq is nybble-encoded according to bam_nt16_table. - */ -typedef struct { - bam1_core_t core; - int l_aux, data_len, m_data; - uint8_t *data; -} bam1_t; - -typedef struct __bam_iter_t *bam_iter_t; - -#define bam1_strand(b) (((b)->core.flag&BAM_FREVERSE) != 0) -#define bam1_mstrand(b) (((b)->core.flag&BAM_FMREVERSE) != 0) - -/*! @function - @abstract Get the CIGAR array - @param b pointer to an alignment - @return pointer to the CIGAR array - - @discussion In the CIGAR array, each element is a 32-bit integer. The - lower 4 bits gives a CIGAR operation and the higher 28 bits keep the - length of a CIGAR. - */ -#define bam1_cigar(b) ((uint32_t*)((b)->data + (b)->core.l_qname)) - -/*! @function - @abstract Get the name of the query - @param b pointer to an alignment - @return pointer to the name string, null terminated - */ -#define bam1_qname(b) ((char*)((b)->data)) - -/*! @function - @abstract Get query sequence - @param b pointer to an alignment - @return pointer to sequence - - @discussion Each base is encoded in 4 bits: 1 for A, 2 for C, 4 for G, - 8 for T and 15 for N. Two bases are packed in one byte with the base - at the higher 4 bits having smaller coordinate on the read. It is - recommended to use bam1_seqi() macro to get the base. - */ -#define bam1_seq(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname) - -/*! @function - @abstract Get query quality - @param b pointer to an alignment - @return pointer to quality string - */ -#define bam1_qual(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (((b)->core.l_qseq + 1)>>1)) - -/*! @function - @abstract Get a base on read - @param s Query sequence returned by bam1_seq() - @param i The i-th position, 0-based - @return 4-bit integer representing the base. - */ -//#define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf) -#define bam1_seqi(s, i) ((s)[(i)>>1] >> ((~(i)&1)<<2) & 0xf) - -#define bam1_seq_seti(s, i, c) ( (s)[(i)>>1] = ((s)[(i)>>1] & 0xf<<(((i)&1)<<2)) | (c)<<((~(i)&1)<<2) ) - -/*! @function - @abstract Get query sequence and quality - @param b pointer to an alignment - @return pointer to the concatenated auxiliary data - */ -#define bam1_aux(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2) - -#ifndef kroundup32 -/*! @function - @abstract Round an integer to the next closest power-2 integer. - @param x integer to be rounded (in place) - @discussion x will be modified. - */ -#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) -#endif - -/*! - @abstract Whether the machine is big-endian; modified only in - bam_header_init(). - */ -extern int bam_is_be; - -/*! - @abstract Verbose level between 0 and 3; 0 is supposed to disable all - debugging information, though this may not have been implemented. - */ -extern int bam_verbose; - -extern int bam_no_B; - -/*! @abstract Table for converting a nucleotide character to the 4-bit encoding. */ -extern unsigned char bam_nt16_table[256]; - -/*! @abstract Table for converting a 4-bit encoded nucleotide to a letter. */ -extern char *bam_nt16_rev_table; - -extern char bam_nt16_nt4_table[]; - -#ifdef __cplusplus -extern "C" { -#endif - - /********************* - * Low-level SAM I/O * - *********************/ - - /*! @abstract TAM file handler */ - typedef struct __tamFile_t *tamFile; - - /*! - @abstract Open a SAM file for reading, either uncompressed or compressed by gzip/zlib. - @param fn SAM file name - @return SAM file handler - */ - tamFile sam_open(const char *fn); - - /*! - @abstract Close a SAM file handler - @param fp SAM file handler - */ - void sam_close(tamFile fp); - - /*! - @abstract Read one alignment from a SAM file handler - @param fp SAM file handler - @param header header information (ordered names of chromosomes) - @param b read alignment; all members in b will be updated - @return 0 if successful; otherwise negative - */ - int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b); - - /*! - @abstract Read header information from a TAB-delimited list file. - @param fn_list file name for the list - @return a pointer to the header structure - - @discussion Each line in this file consists of chromosome name and - the length of chromosome. - */ - bam_header_t *sam_header_read2(const char *fn_list); - - /*! - @abstract Read header from a SAM file (if present) - @param fp SAM file handler - @return pointer to header struct; 0 if no @SQ lines available - */ - bam_header_t *sam_header_read(tamFile fp); - - /*! - @abstract Parse @SQ lines a update a header struct - @param h pointer to the header struct to be updated - @return number of target sequences - - @discussion bam_header_t::{n_targets,target_len,target_name} will - be destroyed in the first place. - */ - int sam_header_parse(bam_header_t *h); - int32_t bam_get_tid(const bam_header_t *header, const char *seq_name); - - /*! - @abstract Parse @RG lines a update a header struct - @param h pointer to the header struct to be updated - @return number of @RG lines - - @discussion bam_header_t::rg2lib will be destroyed in the first - place. - */ - int sam_header_parse_rg(bam_header_t *h); - -#define sam_write1(header, b) bam_view1(header, b) - - - /******************************** - * APIs for string dictionaries * - ********************************/ - - int bam_strmap_put(void *strmap, const char *rg, const char *lib); - const char *bam_strmap_get(const void *strmap, const char *rg); - void *bam_strmap_dup(const void*); - void *bam_strmap_init(); - void bam_strmap_destroy(void *strmap); - - - /********************* - * Low-level BAM I/O * - *********************/ - - /*! - @abstract Initialize a header structure. - @return the pointer to the header structure - - @discussion This function also modifies the global variable - bam_is_be. - */ - bam_header_t *bam_header_init(); - - /*! - @abstract Destroy a header structure. - @param header pointer to the header - */ - void bam_header_destroy(bam_header_t *header); - - /*! - @abstract Read a header structure from BAM. - @param fp BAM file handler, opened by bam_open() - @return pointer to the header structure - - @discussion The file position indicator must be placed at the - beginning of the file. Upon success, the position indicator will - be set at the start of the first alignment. - */ - bam_header_t *bam_header_read(bamFile fp); - - /*! - @abstract Write a header structure to BAM. - @param fp BAM file handler - @param header pointer to the header structure - @return always 0 currently - */ - int bam_header_write(bamFile fp, const bam_header_t *header); - - /*! - @abstract Read an alignment from BAM. - @param fp BAM file handler - @param b read alignment; all members are updated. - @return number of bytes read from the file - - @discussion The file position indicator must be - placed right before an alignment. Upon success, this function - will set the position indicator to the start of the next - alignment. This function is not affected by the machine - endianness. - */ - int bam_read1(bamFile fp, bam1_t *b); - - int bam_remove_B(bam1_t *b); - - /*! - @abstract Write an alignment to BAM. - @param fp BAM file handler - @param c pointer to the bam1_core_t structure - @param data_len total length of variable size data related to - the alignment - @param data pointer to the concatenated data - @return number of bytes written to the file - - @discussion This function is not affected by the machine - endianness. - */ - int bam_write1_core(bamFile fp, const bam1_core_t *c, int data_len, uint8_t *data); - - /*! - @abstract Write an alignment to BAM. - @param fp BAM file handler - @param b alignment to write - @return number of bytes written to the file - - @abstract It is equivalent to: - bam_write1_core(fp, &b->core, b->data_len, b->data) - */ - int bam_write1(bamFile fp, const bam1_t *b); - - /*! @function - @abstract Initiate a pointer to bam1_t struct - */ -#define bam_init1() ((bam1_t*)calloc(1, sizeof(bam1_t))) - - /*! @function - @abstract Free the memory allocated for an alignment. - @param b pointer to an alignment - */ -#define bam_destroy1(b) do { \ - if (b) { free((b)->data); free(b); } \ - } while (0) - - /*! - @abstract Format a BAM record in the SAM format - @param header pointer to the header structure - @param b alignment to print - @return a pointer to the SAM string - */ - char *bam_format1(const bam_header_t *header, const bam1_t *b); - - char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int of); - - /*! - @abstract Check whether a BAM record is plausibly valid - @param header associated header structure, or NULL if unavailable - @param b alignment to validate - @return 0 if the alignment is invalid; non-zero otherwise - - @discussion Simple consistency check of some of the fields of the - alignment record. If the header is provided, several additional checks - are made. Not all fields are checked, so a non-zero result is not a - guarantee that the record is valid. However it is usually good enough - to detect when bam_seek() has been called with a virtual file offset - that is not the offset of an alignment record. - */ - int bam_validate1(const bam_header_t *header, const bam1_t *b); - - const char *bam_get_library(bam_header_t *header, const bam1_t *b); - - - /*************** - * pileup APIs * - ***************/ - - /*! @typedef - @abstract Structure for one alignment covering the pileup position. - @field b pointer to the alignment - @field qpos position of the read base at the pileup site, 0-based - @field indel indel length; 0 for no indel, positive for ins and negative for del - @field is_del 1 iff the base on the padded read is a deletion - @field level the level of the read in the "viewer" mode - - @discussion See also bam_plbuf_push() and bam_lplbuf_push(). The - difference between the two functions is that the former does not - set bam_pileup1_t::level, while the later does. Level helps the - implementation of alignment viewers, but calculating this has some - overhead. - */ - typedef struct { - bam1_t *b; - int32_t qpos; - int indel, level; - uint32_t is_del:1, is_head:1, is_tail:1, is_refskip:1, aux:28; - } bam_pileup1_t; - - typedef int (*bam_plp_auto_f)(void *data, bam1_t *b); - - struct __bam_plp_t; - typedef struct __bam_plp_t *bam_plp_t; - - bam_plp_t bam_plp_init(bam_plp_auto_f func, void *data); - int bam_plp_push(bam_plp_t iter, const bam1_t *b); - const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp); - const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp); - void bam_plp_set_mask(bam_plp_t iter, int mask); - void bam_plp_set_maxcnt(bam_plp_t iter, int maxcnt); - void bam_plp_reset(bam_plp_t iter); - void bam_plp_destroy(bam_plp_t iter); - - struct __bam_mplp_t; - typedef struct __bam_mplp_t *bam_mplp_t; - - bam_mplp_t bam_mplp_init(int n, bam_plp_auto_f func, void **data); - void bam_mplp_destroy(bam_mplp_t iter); - void bam_mplp_set_maxcnt(bam_mplp_t iter, int maxcnt); - int bam_mplp_auto(bam_mplp_t iter, int *_tid, int *_pos, int *n_plp, const bam_pileup1_t **plp); - - /*! @typedef - @abstract Type of function to be called by bam_plbuf_push(). - @param tid chromosome ID as is defined in the header - @param pos start coordinate of the alignment, 0-based - @param n number of elements in pl array - @param pl array of alignments - @param data user provided data - @discussion See also bam_plbuf_push(), bam_plbuf_init() and bam_pileup1_t. - */ - typedef int (*bam_pileup_f)(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data); - - typedef struct { - bam_plp_t iter; - bam_pileup_f func; - void *data; - } bam_plbuf_t; - - void bam_plbuf_set_mask(bam_plbuf_t *buf, int mask); - void bam_plbuf_reset(bam_plbuf_t *buf); - bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data); - void bam_plbuf_destroy(bam_plbuf_t *buf); - int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf); - - int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data); - - struct __bam_lplbuf_t; - typedef struct __bam_lplbuf_t bam_lplbuf_t; - - void bam_lplbuf_reset(bam_lplbuf_t *buf); - - /*! @abstract bam_plbuf_init() equivalent with level calculated. */ - bam_lplbuf_t *bam_lplbuf_init(bam_pileup_f func, void *data); - - /*! @abstract bam_plbuf_destroy() equivalent with level calculated. */ - void bam_lplbuf_destroy(bam_lplbuf_t *tv); - - /*! @abstract bam_plbuf_push() equivalent with level calculated. */ - int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *buf); - - - /********************* - * BAM indexing APIs * - *********************/ - - struct __bam_index_t; - typedef struct __bam_index_t bam_index_t; - - /*! - @abstract Build index for a BAM file. - @discussion Index file "fn.bai" will be created. - @param fn name of the BAM file - @return always 0 currently - */ - int bam_index_build(const char *fn); - - /*! - @abstract Load index from file "fn.bai". - @param fn name of the BAM file (NOT the index file) - @return pointer to the index structure - */ - bam_index_t *bam_index_load(const char *fn); - - /*! - @abstract Destroy an index structure. - @param idx pointer to the index structure - */ - void bam_index_destroy(bam_index_t *idx); - - /*! @typedef - @abstract Type of function to be called by bam_fetch(). - @param b the alignment - @param data user provided data - */ - typedef int (*bam_fetch_f)(const bam1_t *b, void *data); - - /*! - @abstract Retrieve the alignments that are overlapped with the - specified region. - - @discussion A user defined function will be called for each - retrieved alignment ordered by its start position. - - @param fp BAM file handler - @param idx pointer to the alignment index - @param tid chromosome ID as is defined in the header - @param beg start coordinate, 0-based - @param end end coordinate, 0-based - @param data user provided data (will be transferred to func) - @param func user defined function - */ - int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func); - - bam_iter_t bam_iter_query(const bam_index_t *idx, int tid, int beg, int end); - int bam_iter_read(bamFile fp, bam_iter_t iter, bam1_t *b); - void bam_iter_destroy(bam_iter_t iter); - - /*! - @abstract Parse a region in the format: "chr2:100,000-200,000". - @discussion bam_header_t::hash will be initialized if empty. - @param header pointer to the header structure - @param str string to be parsed - @param ref_id the returned chromosome ID - @param begin the returned start coordinate - @param end the returned end coordinate - @return 0 on success; -1 on failure - */ - int bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *begin, int *end); - - - /************************** - * APIs for optional tags * - **************************/ - - /*! - @abstract Retrieve data of a tag - @param b pointer to an alignment struct - @param tag two-character tag to be retrieved - - @return pointer to the type and data. The first character is the - type that can be 'iIsScCdfAZH'. - - @discussion Use bam_aux2?() series to convert the returned data to - the corresponding type. - */ - uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]); - - int32_t bam_aux2i(const uint8_t *s); - float bam_aux2f(const uint8_t *s); - double bam_aux2d(const uint8_t *s); - char bam_aux2A(const uint8_t *s); - char *bam_aux2Z(const uint8_t *s); - - int bam_aux_del(bam1_t *b, uint8_t *s); - void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data); - uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]); // an alias of bam_aux_get() - - - /***************** - * Miscellaneous * - *****************/ - - /*! - @abstract Calculate the rightmost coordinate of an alignment on the - reference genome. - - @param c pointer to the bam1_core_t structure - @param cigar the corresponding CIGAR array (from bam1_t::cigar) - @return the rightmost coordinate, 0-based - */ - uint32_t bam_calend(const bam1_core_t *c, const uint32_t *cigar); - - /*! - @abstract Calculate the length of the query sequence from CIGAR. - @param c pointer to the bam1_core_t structure - @param cigar the corresponding CIGAR array (from bam1_t::cigar) - @return length of the query sequence - */ - int32_t bam_cigar2qlen(const bam1_core_t *c, const uint32_t *cigar); - -#ifdef __cplusplus -} -#endif - -/*! - @abstract Calculate the minimum bin that contains a region [beg,end). - @param beg start of the region, 0-based - @param end end of the region, 0-based - @return bin - */ -static inline int bam_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; -} - -/*! - @abstract Copy an alignment - @param bdst destination alignment struct - @param bsrc source alignment struct - @return pointer to the destination alignment struct - */ -static inline bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc) -{ - uint8_t *data = bdst->data; - int m_data = bdst->m_data; // backup data and m_data - if (m_data < bsrc->data_len) { // double the capacity - m_data = bsrc->data_len; kroundup32(m_data); - data = (uint8_t*)realloc(data, m_data); - } - memcpy(data, bsrc->data, bsrc->data_len); // copy var-len data - *bdst = *bsrc; // copy the rest - // restore the backup - bdst->m_data = m_data; - bdst->data = data; - return bdst; -} - -/*! - @abstract Duplicate an alignment - @param src source alignment struct - @return pointer to the destination alignment struct - */ -static inline bam1_t *bam_dup1(const bam1_t *src) -{ - bam1_t *b; - b = bam_init1(); - *b = *src; - b->m_data = b->data_len; - b->data = (uint8_t*)calloc(b->data_len, 1); - memcpy(b->data, src->data, b->data_len); - return b; -} - -static inline int bam_aux_type2size(int x) -{ - if (x == 'C' || x == 'c' || x == 'A') return 1; - else if (x == 'S' || x == 's') return 2; - else if (x == 'I' || x == 'i' || x == 'f' || x == 'F') return 4; - else return 0; -} - -/********************************* - *** Compatibility with htslib *** - *********************************/ - -typedef bam_header_t bam_hdr_t; - -#define bam_get_qname(b) bam1_qname(b) -#define bam_get_cigar(b) bam1_cigar(b) - -#define bam_hdr_read(fp) bam_header_read(fp) -#define bam_hdr_write(fp, h) bam_header_write(fp, h) -#define bam_hdr_destroy(fp) bam_header_destroy(fp) - -#endif diff --git a/samlib/bam_aux.c b/samlib/bam_aux.c deleted file mode 100644 index 4bbf975..0000000 --- a/samlib/bam_aux.c +++ /dev/null @@ -1,217 +0,0 @@ -#include -#include "bam.h" -#include "khash.h" -typedef char *str_p; -KHASH_MAP_INIT_STR(s, int) -KHASH_MAP_INIT_STR(r2l, str_p) - -void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data) -{ - int ori_len = b->data_len; - b->data_len += 3 + len; - b->l_aux += 3 + len; - if (b->m_data < b->data_len) { - b->m_data = b->data_len; - kroundup32(b->m_data); - b->data = (uint8_t*)realloc(b->data, b->m_data); - } - b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1]; - b->data[ori_len + 2] = type; - memcpy(b->data + ori_len + 3, data, len); -} - -uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]) -{ - return bam_aux_get(b, tag); -} - -#define __skip_tag(s) do { \ - int type = toupper(*(s)); \ - ++(s); \ - if (type == 'Z' || type == 'H') { while (*(s)) ++(s); ++(s); } \ - else if (type == 'B') (s) += 5 + bam_aux_type2size(*(s)) * (*(int32_t*)((s)+1)); \ - else (s) += bam_aux_type2size(type); \ - } while(0) - -uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]) -{ - uint8_t *s; - int y = tag[0]<<8 | tag[1]; - s = bam1_aux(b); - while (s < b->data + b->data_len) { - int x = (int)s[0]<<8 | s[1]; - s += 2; - if (x == y) return s; - __skip_tag(s); - } - return 0; -} -// s MUST BE returned by bam_aux_get() -int bam_aux_del(bam1_t *b, uint8_t *s) -{ - uint8_t *p, *aux; - aux = bam1_aux(b); - p = s - 2; - __skip_tag(s); - memmove(p, s, b->l_aux - (s - aux)); - b->data_len -= s - p; - b->l_aux -= s - p; - return 0; -} - -int bam_aux_drop_other(bam1_t *b, uint8_t *s) -{ - if (s) { - uint8_t *p, *aux; - aux = bam1_aux(b); - p = s - 2; - __skip_tag(s); - memmove(aux, p, s - p); - b->data_len -= b->l_aux - (s - p); - b->l_aux = s - p; - } else { - b->data_len -= b->l_aux; - b->l_aux = 0; - } - return 0; -} - -void bam_init_header_hash(bam_header_t *header) -{ - if (header->hash == 0) { - int ret, i; - khiter_t iter; - khash_t(s) *h; - header->hash = h = kh_init(s); - for (i = 0; i < header->n_targets; ++i) { - iter = kh_put(s, h, header->target_name[i], &ret); - kh_value(h, iter) = i; - } - } -} - -void bam_destroy_header_hash(bam_header_t *header) -{ - if (header->hash) - kh_destroy(s, (khash_t(s)*)header->hash); -} - -int32_t bam_get_tid(const bam_header_t *header, const char *seq_name) -{ - khint_t k; - khash_t(s) *h = (khash_t(s)*)header->hash; - k = kh_get(s, h, seq_name); - return k == kh_end(h)? -1 : kh_value(h, k); -} - -int bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *beg, int *end) -{ - char *s; - int i, l, k, name_end; - khiter_t iter; - khash_t(s) *h; - - bam_init_header_hash(header); - h = (khash_t(s)*)header->hash; - - *ref_id = *beg = *end = -1; - name_end = l = strlen(str); - s = (char*)malloc(l+1); - // remove space - for (i = k = 0; i < l; ++i) - if (!isspace(str[i])) s[k++] = str[i]; - s[k] = 0; l = k; - // determine the sequence name - for (i = l - 1; i >= 0; --i) if (s[i] == ':') break; // look for colon from the end - if (i >= 0) name_end = i; - if (name_end < l) { // check if this is really the end - int n_hyphen = 0; - for (i = name_end + 1; i < l; ++i) { - if (s[i] == '-') ++n_hyphen; - else if (!isdigit(s[i]) && s[i] != ',') break; - } - if (i < l || n_hyphen > 1) name_end = l; // malformated region string; then take str as the name - s[name_end] = 0; - iter = kh_get(s, h, s); - if (iter == kh_end(h)) { // cannot find the sequence name - iter = kh_get(s, h, str); // try str as the name - if (iter == kh_end(h)) { - if (bam_verbose >= 2) fprintf(stderr, "[%s] fail to determine the sequence name.\n", __func__); - free(s); return -1; - } else s[name_end] = ':', name_end = l; - } - } else iter = kh_get(s, h, str); - if (iter == kh_end(h)) { - free(s); - return -1; - } - *ref_id = kh_val(h, iter); - // parse the interval - if (name_end < l) { - for (i = k = name_end + 1; i < l; ++i) - if (s[i] != ',') s[k++] = s[i]; - s[k] = 0; - *beg = atoi(s + name_end + 1); - for (i = name_end + 1; i != k; ++i) if (s[i] == '-') break; - *end = i < k? atoi(s + i + 1) : 1<<29; - if (*beg > 0) --*beg; - } else *beg = 0, *end = 1<<29; - free(s); - return *beg <= *end? 0 : -1; -} - -int32_t bam_aux2i(const uint8_t *s) -{ - int type; - if (s == 0) return 0; - type = *s++; - if (type == 'c') return (int32_t)*(int8_t*)s; - else if (type == 'C') return (int32_t)*(uint8_t*)s; - else if (type == 's') return (int32_t)*(int16_t*)s; - else if (type == 'S') return (int32_t)*(uint16_t*)s; - else if (type == 'i' || type == 'I') return *(int32_t*)s; - else return 0; -} - -float bam_aux2f(const uint8_t *s) -{ - int type; - type = *s++; - if (s == 0) return 0.0; - if (type == 'f') return *(float*)s; - else return 0.0; -} - -double bam_aux2d(const uint8_t *s) -{ - int type; - type = *s++; - if (s == 0) return 0.0; - if (type == 'd') return *(double*)s; - else return 0.0; -} - -char bam_aux2A(const uint8_t *s) -{ - int type; - type = *s++; - if (s == 0) return 0; - if (type == 'A') return *(char*)s; - else return 0; -} - -char *bam_aux2Z(const uint8_t *s) -{ - int type; - type = *s++; - if (s == 0) return 0; - if (type == 'Z' || type == 'H') return (char*)s; - else return 0; -} - -#ifdef _WIN32 -double drand48() -{ - return (double)rand() / RAND_MAX; -} -#endif diff --git a/samlib/bam_endian.h b/samlib/bam_endian.h deleted file mode 100644 index 0fc74a8..0000000 --- a/samlib/bam_endian.h +++ /dev/null @@ -1,42 +0,0 @@ -#ifndef BAM_ENDIAN_H -#define BAM_ENDIAN_H - -#include - -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/samlib/bam_import.c b/samlib/bam_import.c deleted file mode 100644 index da2bf94..0000000 --- a/samlib/bam_import.c +++ /dev/null @@ -1,489 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#ifdef _WIN32 -#include -#endif -#include "kstring.h" -#include "bam.h" -#include "sam_header.h" -#include "kseq.h" -#include "khash.h" - -KSTREAM_INIT(gzFile, gzread, 16384) -KHASH_MAP_INIT_STR(ref, uint64_t) - -void bam_init_header_hash(bam_header_t *header); -void bam_destroy_header_hash(bam_header_t *header); -int32_t bam_get_tid(const bam_header_t *header, const char *seq_name); - -unsigned char bam_nt16_table[256] = { - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 1, 2, 4, 8, 15,15,15,15, 15,15,15,15, 15, 0 /*=*/,15,15, - 15, 1,14, 2, 13,15,15, 4, 11,15,15,12, 15, 3,15,15, - 15,15, 5, 6, 8,15, 7, 9, 15,10,15,15, 15,15,15,15, - 15, 1,14, 2, 13,15,15, 4, 11,15,15,12, 15, 3,15,15, - 15,15, 5, 6, 8,15, 7, 9, 15,10,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, - 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15 -}; - -unsigned short bam_char2flag_table[256] = { - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,BAM_FREAD1,BAM_FREAD2,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - BAM_FPROPER_PAIR,0,BAM_FMREVERSE,0, 0,BAM_FMUNMAP,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, BAM_FDUP,0,BAM_FQCFAIL,0, 0,0,0,0, 0,0,0,0, - BAM_FPAIRED,0,BAM_FREVERSE,BAM_FSECONDARY, 0,BAM_FUNMAP,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, - 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 -}; - -char *bam_nt16_rev_table = "=ACMGRSVTWYHKDBN"; - -struct __tamFile_t { - gzFile fp; - kstream_t *ks; - kstring_t *str; - uint64_t n_lines; - int is_first; -}; - -char **__bam_get_lines(const char *fn, int *_n) // for bam_plcmd.c only -{ - char **list = 0, *s; - int n = 0, dret, m = 0; - gzFile fp = (strcmp(fn, "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(fn, "r"); - kstream_t *ks; - kstring_t *str; - str = (kstring_t*)calloc(1, sizeof(kstring_t)); - ks = ks_init(fp); - while (ks_getuntil(ks, '\n', str, &dret) > 0) { - if (n == m) { - m = m? m << 1 : 16; - list = (char**)realloc(list, m * sizeof(char*)); - } - if (str->s[str->l-1] == '\r') - str->s[--str->l] = '\0'; - s = list[n++] = (char*)calloc(str->l + 1, 1); - strcpy(s, str->s); - } - ks_destroy(ks); - gzclose(fp); - free(str->s); free(str); - *_n = n; - return list; -} - -static bam_header_t *hash2header(const kh_ref_t *hash) -{ - bam_header_t *header; - khiter_t k; - header = bam_header_init(); - header->n_targets = kh_size(hash); - header->target_name = (char**)calloc(kh_size(hash), sizeof(char*)); - header->target_len = (uint32_t*)calloc(kh_size(hash), 4); - for (k = kh_begin(hash); k != kh_end(hash); ++k) { - if (kh_exist(hash, k)) { - int i = (int)kh_value(hash, k); - header->target_name[i] = (char*)kh_key(hash, k); - header->target_len[i] = kh_value(hash, k)>>32; - } - } - bam_init_header_hash(header); - return header; -} -bam_header_t *sam_header_read2(const char *fn) -{ - bam_header_t *header; - int c, dret, ret, error = 0; - gzFile fp; - kstream_t *ks; - kstring_t *str; - kh_ref_t *hash; - khiter_t k; - if (fn == 0) return 0; - fp = (strcmp(fn, "-") == 0)? gzdopen(fileno(stdin), "r") : gzopen(fn, "r"); - if (fp == 0) return 0; - hash = kh_init(ref); - ks = ks_init(fp); - str = (kstring_t*)calloc(1, sizeof(kstring_t)); - while (ks_getuntil(ks, 0, str, &dret) > 0) { - char *s = strdup(str->s); - int len, i; - i = kh_size(hash); - ks_getuntil(ks, 0, str, &dret); - len = atoi(str->s); - k = kh_put(ref, hash, s, &ret); - if (ret == 0) { - fprintf(stderr, "[sam_header_read2] duplicated sequence name: %s\n", s); - error = 1; - } - kh_value(hash, k) = (uint64_t)len<<32 | i; - if (dret != '\n') - while ((c = ks_getc(ks)) != '\n' && c != -1); - } - ks_destroy(ks); - gzclose(fp); - free(str->s); free(str); - fprintf(stderr, "[sam_header_read2] %d sequences loaded.\n", kh_size(hash)); - if (error) return 0; - header = hash2header(hash); - kh_destroy(ref, hash); - return header; -} -static inline uint8_t *alloc_data(bam1_t *b, int size) -{ - if (b->m_data < size) { - b->m_data = size; - kroundup32(b->m_data); - b->data = (uint8_t*)realloc(b->data, b->m_data); - } - return b->data; -} -static inline void parse_error(int64_t n_lines, const char * __restrict msg) -{ - fprintf(stderr, "Parse error at line %lld: %s\n", (long long)n_lines, msg); - abort(); -} -static inline void append_text(bam_header_t *header, kstring_t *str) -{ - size_t x = header->l_text, y = header->l_text + str->l + 2; // 2 = 1 byte dret + 1 byte null - kroundup32(x); kroundup32(y); - if (x < y) - { - header->n_text = y; - header->text = (char*)realloc(header->text, y); - if ( !header->text ) - { - fprintf(stderr,"realloc failed to alloc %ld bytes\n", y); - abort(); - } - } - // Sanity check - if ( header->l_text+str->l+1 >= header->n_text ) - { - fprintf(stderr,"append_text FIXME: %ld>=%ld, x=%ld,y=%ld\n", header->l_text+str->l+1,(long)header->n_text,x,y); - abort(); - } - strncpy(header->text + header->l_text, str->s, str->l+1); // we cannot use strcpy() here. - header->l_text += str->l + 1; - header->text[header->l_text] = 0; -} - -int sam_header_parse(bam_header_t *h) -{ - char **tmp; - int i; - free(h->target_len); free(h->target_name); - h->n_targets = 0; h->target_len = 0; h->target_name = 0; - if (h->l_text < 3) return 0; - if (h->dict == 0) h->dict = sam_header_parse2(h->text); - tmp = sam_header2list(h->dict, "SQ", "SN", &h->n_targets); - if (h->n_targets == 0) return 0; - h->target_name = calloc(h->n_targets, sizeof(void*)); - for (i = 0; i < h->n_targets; ++i) - h->target_name[i] = strdup(tmp[i]); - free(tmp); - tmp = sam_header2list(h->dict, "SQ", "LN", &h->n_targets); - h->target_len = calloc(h->n_targets, 4); - for (i = 0; i < h->n_targets; ++i) - h->target_len[i] = atoi(tmp[i]); - free(tmp); - return h->n_targets; -} - -bam_header_t *sam_header_read(tamFile fp) -{ - int ret, dret; - bam_header_t *header = bam_header_init(); - kstring_t *str = fp->str; - while ((ret = ks_getuntil(fp->ks, KS_SEP_TAB, str, &dret)) >= 0 && str->s[0] == '@') { // skip header - str->s[str->l] = dret; // note that str->s is NOT null terminated!! - append_text(header, str); - if (dret != '\n') { - ret = ks_getuntil(fp->ks, '\n', str, &dret); - str->s[str->l] = '\n'; // NOT null terminated!! - append_text(header, str); - } - ++fp->n_lines; - } - sam_header_parse(header); - bam_init_header_hash(header); - fp->is_first = 1; - return header; -} - -int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b) -{ - int ret, doff, doff0, dret, z = 0; - bam1_core_t *c = &b->core; - kstring_t *str = fp->str; - kstream_t *ks = fp->ks; - - if (fp->is_first) { - fp->is_first = 0; - ret = str->l; - } else { - do { // special consideration for empty lines - ret = ks_getuntil(fp->ks, KS_SEP_TAB, str, &dret); - if (ret >= 0) z += str->l + 1; - } while (ret == 0); - } - if (ret < 0) return -1; - ++fp->n_lines; - doff = 0; - - { // name - c->l_qname = strlen(str->s) + 1; - memcpy(alloc_data(b, doff + c->l_qname) + doff, str->s, c->l_qname); - doff += c->l_qname; - } - { // flag - long flag; - char *s; - ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; - flag = strtol((char*)str->s, &s, 0); - if (*s) { // not the end of the string - flag = 0; - for (s = str->s; *s; ++s) - flag |= bam_char2flag_table[(int)*s]; - } - c->flag = flag; - } - { // tid, pos, qual - ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; c->tid = bam_get_tid(header, str->s); - if (c->tid < 0 && strcmp(str->s, "*")) { - if (header->n_targets == 0) { - fprintf(stderr, "[sam_read1] missing header? Abort!\n"); - exit(1); - } else fprintf(stderr, "[sam_read1] reference '%s' is recognized as '*'.\n", str->s); - } - ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; c->pos = isdigit(str->s[0])? atoi(str->s) - 1 : -1; - ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; c->qual = isdigit(str->s[0])? atoi(str->s) : 0; - if (ret < 0) return -2; - } - { // cigar - char *s, *t; - int i, op; - long x; - c->n_cigar = 0; - if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -3; - z += str->l + 1; - if (str->s[0] != '*') { - uint32_t *cigar; - for (s = str->s; *s; ++s) { - if ((isalpha(*s)) || (*s=='=')) ++c->n_cigar; - else if (!isdigit(*s)) parse_error(fp->n_lines, "invalid CIGAR character"); - } - b->data = alloc_data(b, doff + c->n_cigar * 4); - cigar = bam1_cigar(b); - for (i = 0, s = str->s; i != c->n_cigar; ++i) { - x = strtol(s, &t, 10); - op = toupper(*t); - if (op == 'M') op = BAM_CMATCH; - else if (op == 'I') op = BAM_CINS; - else if (op == 'D') op = BAM_CDEL; - else if (op == 'N') op = BAM_CREF_SKIP; - else if (op == 'S') op = BAM_CSOFT_CLIP; - else if (op == 'H') op = BAM_CHARD_CLIP; - else if (op == 'P') op = BAM_CPAD; - else if (op == '=') op = BAM_CEQUAL; - else if (op == 'X') op = BAM_CDIFF; - else if (op == 'B') op = BAM_CBACK; - else parse_error(fp->n_lines, "invalid CIGAR operation"); - s = t + 1; - cigar[i] = bam_cigar_gen(x, op); - } - if (*s) parse_error(fp->n_lines, "unmatched CIGAR operation"); - c->bin = bam_reg2bin(c->pos, bam_calend(c, cigar)); - doff += c->n_cigar * 4; - } else { - if (!(c->flag&BAM_FUNMAP)) { - fprintf(stderr, "Parse warning at line %lld: mapped sequence without CIGAR\n", (long long)fp->n_lines); - c->flag |= BAM_FUNMAP; - } - c->bin = bam_reg2bin(c->pos, c->pos + 1); - } - } - { // mtid, mpos, isize - ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; - c->mtid = strcmp(str->s, "=")? bam_get_tid(header, str->s) : c->tid; - ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; - c->mpos = isdigit(str->s[0])? atoi(str->s) - 1 : -1; - ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; - c->isize = (str->s[0] == '-' || isdigit(str->s[0]))? atoi(str->s) : 0; - if (ret < 0) return -4; - } - { // seq and qual - int i; - uint8_t *p = 0; - if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -5; // seq - z += str->l + 1; - if (strcmp(str->s, "*")) { - c->l_qseq = strlen(str->s); - if (c->n_cigar && c->l_qseq != (int32_t)bam_cigar2qlen(c, bam1_cigar(b))) { - fprintf(stderr, "Line %ld, sequence length %i vs %i from CIGAR\n", - (long)fp->n_lines, c->l_qseq, (int32_t)bam_cigar2qlen(c, bam1_cigar(b))); - parse_error(fp->n_lines, "CIGAR and sequence length are inconsistent"); - } - p = (uint8_t*)alloc_data(b, doff + c->l_qseq + (c->l_qseq+1)/2) + doff; - memset(p, 0, (c->l_qseq+1)/2); - for (i = 0; i < c->l_qseq; ++i) - p[i/2] |= bam_nt16_table[(int)str->s[i]] << 4*(1-i%2); - } else c->l_qseq = 0; - if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -6; // qual - z += str->l + 1; - if (strcmp(str->s, "*") && c->l_qseq != strlen(str->s)) - parse_error(fp->n_lines, "sequence and quality are inconsistent"); - p += (c->l_qseq+1)/2; - if (strcmp(str->s, "*") == 0) for (i = 0; i < c->l_qseq; ++i) p[i] = 0xff; - else for (i = 0; i < c->l_qseq; ++i) p[i] = str->s[i] - 33; - doff += c->l_qseq + (c->l_qseq+1)/2; - } - doff0 = doff; - if (dret != '\n' && dret != '\r') { // aux - while (ks_getuntil(ks, KS_SEP_TAB, str, &dret) >= 0) { - uint8_t *s, type, key[2]; - z += str->l + 1; - if (str->l < 6 || str->s[2] != ':' || str->s[4] != ':') - parse_error(fp->n_lines, "missing colon in auxiliary data"); - key[0] = str->s[0]; key[1] = str->s[1]; - type = str->s[3]; - s = alloc_data(b, doff + 3) + doff; - s[0] = key[0]; s[1] = key[1]; s += 2; doff += 2; - if (type == 'A' || type == 'a' || type == 'c' || type == 'C') { // c and C for backward compatibility - s = alloc_data(b, doff + 2) + doff; - *s++ = 'A'; *s = str->s[5]; - doff += 2; - } else if (type == 'I' || type == 'i') { - long long x; - s = alloc_data(b, doff + 5) + doff; - x = (long long)atoll(str->s + 5); - if (x < 0) { - if (x >= -127) { - *s++ = 'c'; *(int8_t*)s = (int8_t)x; - s += 1; doff += 2; - } else if (x >= -32767) { - *s++ = 's'; *(int16_t*)s = (int16_t)x; - s += 2; doff += 3; - } else { - *s++ = 'i'; *(int32_t*)s = (int32_t)x; - s += 4; doff += 5; - if (x < -2147483648ll) - fprintf(stderr, "Parse warning at line %lld: integer %lld is out of range.", - (long long)fp->n_lines, x); - } - } else { - if (x <= 255) { - *s++ = 'C'; *s++ = (uint8_t)x; - doff += 2; - } else if (x <= 65535) { - *s++ = 'S'; *(uint16_t*)s = (uint16_t)x; - s += 2; doff += 3; - } else { - *s++ = 'I'; *(uint32_t*)s = (uint32_t)x; - s += 4; doff += 5; - if (x > 4294967295ll) - fprintf(stderr, "Parse warning at line %lld: integer %lld is out of range.", - (long long)fp->n_lines, x); - } - } - } else if (type == 'f') { - s = alloc_data(b, doff + 5) + doff; - *s++ = 'f'; - *(float*)s = (float)atof(str->s + 5); - s += 4; doff += 5; - } else if (type == 'd') { - s = alloc_data(b, doff + 9) + doff; - *s++ = 'd'; - *(float*)s = (float)atof(str->s + 9); - s += 8; doff += 9; - } else if (type == 'Z' || type == 'H') { - int size = 1 + (str->l - 5) + 1; - if (type == 'H') { // check whether the hex string is valid - int i; - if ((str->l - 5) % 2 == 1) parse_error(fp->n_lines, "length of the hex string not even"); - for (i = 0; i < str->l - 5; ++i) { - int c = toupper(str->s[5 + i]); - if (!((c >= '0' && c <= '9') || (c >= 'A' && c <= 'F'))) - parse_error(fp->n_lines, "invalid hex character"); - } - } - s = alloc_data(b, doff + size) + doff; - *s++ = type; - memcpy(s, str->s + 5, str->l - 5); - s[str->l - 5] = 0; - doff += size; - } else if (type == 'B') { - int32_t n = 0, Bsize, k = 0, size; - char *p; - if (str->l < 8) parse_error(fp->n_lines, "too few values in aux type B"); - Bsize = bam_aux_type2size(str->s[5]); // the size of each element - for (p = (char*)str->s + 6; *p; ++p) // count the number of elements in the array - if (*p == ',') ++n; - p = str->s + 7; // now p points to the first number in the array - size = 6 + Bsize * n; // total number of bytes allocated to this tag - s = alloc_data(b, doff + 6 * Bsize * n) + doff; // allocate memory - *s++ = 'B'; *s++ = str->s[5]; - memcpy(s, &n, 4); s += 4; // write the number of elements - if (str->s[5] == 'c') while (p < str->s + str->l) ((int8_t*)s)[k++] = (int8_t)strtol(p, &p, 0), ++p; - else if (str->s[5] == 'C') while (p < str->s + str->l) ((uint8_t*)s)[k++] = (uint8_t)strtol(p, &p, 0), ++p; - else if (str->s[5] == 's') while (p < str->s + str->l) ((int16_t*)s)[k++] = (int16_t)strtol(p, &p, 0), ++p; // FIXME: avoid unaligned memory - else if (str->s[5] == 'S') while (p < str->s + str->l) ((uint16_t*)s)[k++] = (uint16_t)strtol(p, &p, 0), ++p; - else if (str->s[5] == 'i') while (p < str->s + str->l) ((int32_t*)s)[k++] = (int32_t)strtol(p, &p, 0), ++p; - else if (str->s[5] == 'I') while (p < str->s + str->l) ((uint32_t*)s)[k++] = (uint32_t)strtol(p, &p, 0), ++p; - else if (str->s[5] == 'f') while (p < str->s + str->l) ((float*)s)[k++] = (float)strtod(p, &p), ++p; - else parse_error(fp->n_lines, "unrecognized array type"); - s += Bsize * n; doff += size; - } else parse_error(fp->n_lines, "unrecognized type"); - if (dret == '\n' || dret == '\r') break; - } - } - b->l_aux = doff - doff0; - b->data_len = doff; - if (bam_no_B) bam_remove_B(b); - return z; -} - -tamFile sam_open(const char *fn) -{ - tamFile fp; - gzFile gzfp = (strcmp(fn, "-") == 0)? gzdopen(fileno(stdin), "rb") : gzopen(fn, "rb"); - if (gzfp == 0) return 0; - fp = (tamFile)calloc(1, sizeof(struct __tamFile_t)); - fp->str = (kstring_t*)calloc(1, sizeof(kstring_t)); - fp->fp = gzfp; - fp->ks = ks_init(fp->fp); - return fp; -} - -void sam_close(tamFile fp) -{ - if (fp) { - ks_destroy(fp->ks); - gzclose(fp->fp); - free(fp->str->s); free(fp->str); - free(fp); - } -} diff --git a/samlib/bam_index.c b/samlib/bam_index.c deleted file mode 100644 index f916e04..0000000 --- a/samlib/bam_index.c +++ /dev/null @@ -1,726 +0,0 @@ -#include -#include -#include "bam.h" -#include "khash.h" -#include "ksort.h" -#include "bam_endian.h" -#ifdef _USE_KNETFILE -#include "knetfile.h" -#endif - -/*! - @header - - Alignment indexing. Before indexing, BAM must be sorted based on the - leftmost coordinate of alignments. In indexing, BAM uses two indices: - a UCSC binning index and a simple linear index. The binning index is - efficient for alignments spanning long distance, while the auxiliary - linear index helps to reduce unnecessary seek calls especially for - short alignments. - - The UCSC binning scheme was suggested by Richard Durbin and Lincoln - Stein and is explained by Kent et al. (2002). In this scheme, each bin - represents a contiguous genomic region which can be fully contained in - another bin; each alignment is associated with a bin which represents - the smallest region containing the entire alignment. The binning - scheme is essentially another representation of R-tree. A distinct bin - uniquely corresponds to a distinct internal node in a R-tree. Bin A is - a child of Bin B if region A is contained in B. - - In BAM, each bin may span 2^29, 2^26, 2^23, 2^20, 2^17 or 2^14 bp. Bin - 0 spans a 512Mbp region, bins 1-8 span 64Mbp, 9-72 8Mbp, 73-584 1Mbp, - 585-4680 128Kbp and bins 4681-37449 span 16Kbp regions. If we want to - find the alignments overlapped with a region [rbeg,rend), we need to - calculate the list of bins that may be overlapped the region and test - the alignments in the bins to confirm the overlaps. If the specified - region is short, typically only a few alignments in six bins need to - be retrieved. The overlapping alignments can be quickly fetched. - - */ - -#define BAM_MIN_CHUNK_GAP 32768 -// 1<<14 is the size of minimum bin. -#define BAM_LIDX_SHIFT 14 - -#define BAM_MAX_BIN 37450 // =(8^6-1)/7+1 - -typedef struct { - uint64_t u, v; -} pair64_t; - -#define pair64_lt(a,b) ((a).u < (b).u) -KSORT_INIT(off, pair64_t, pair64_lt) - -typedef struct { - uint32_t m, n; - pair64_t *list; -} bam_binlist_t; - -typedef struct { - int32_t n, m; - uint64_t *offset; -} bam_lidx_t; - -KHASH_MAP_INIT_INT(i, bam_binlist_t) - -struct __bam_index_t { - int32_t n; - uint64_t n_no_coor; // unmapped reads without coordinate - khash_t(i) **index; - bam_lidx_t *index2; -}; - -// requirement: len <= LEN_MASK -static inline void insert_offset(khash_t(i) *h, int bin, uint64_t beg, uint64_t end) -{ - khint_t k; - bam_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 void insert_offset2(bam_lidx_t *index2, bam1_t *b, uint64_t offset) -{ - int i, beg, end; - beg = b->core.pos >> BAM_LIDX_SHIFT; - end = (bam_calend(&b->core, bam1_cigar(b)) - 1) >> BAM_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; - } - index2->n = end + 1; -} - -static void merge_chunks(bam_index_t *idx) -{ -#if defined(BAM_TRUE_OFFSET) || defined(BAM_VIRTUAL_OFFSET16) - 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) { - bam_binlist_t *p; - if (!kh_exist(index, k) || kh_key(index, k) == BAM_MAX_BIN) continue; - p = &kh_value(index, k); - m = 0; - for (l = 1; l < p->n; ++l) { -#ifdef BAM_TRUE_OFFSET - if (p->list[m].v + BAM_MIN_CHUNK_GAP > p->list[l].u) p->list[m].v = p->list[l].v; -#else - if (p->list[m].v>>16 == p->list[l].u>>16) p->list[m].v = p->list[l].v; -#endif - else p->list[++m] = p->list[l]; - } // ~for(l) - p->n = m + 1; - } // ~for(k) - } // ~for(i) -#endif // defined(BAM_TRUE_OFFSET) || defined(BAM_BGZF) -} - -static void fill_missing(bam_index_t *idx) -{ - int i, j; - for (i = 0; i < idx->n; ++i) { - bam_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]; - } -} - -bam_index_t *bam_index_core(bamFile fp) -{ - bam1_t *b; - bam_header_t *h; - int i, ret; - bam_index_t *idx; - uint32_t last_bin, save_bin; - int32_t last_coor, last_tid, save_tid; - bam1_core_t *c; - uint64_t save_off, last_off, n_mapped, n_unmapped, off_beg, off_end, n_no_coor; - - h = bam_header_read(fp); - if(h == 0) { - fprintf(stderr, "[bam_index_core] Invalid BAM header."); - return NULL; - } - - idx = (bam_index_t*)calloc(1, sizeof(bam_index_t)); - b = (bam1_t*)calloc(1, sizeof(bam1_t)); - c = &b->core; - - idx->n = h->n_targets; - bam_header_destroy(h); - idx->index = (khash_t(i)**)calloc(idx->n, sizeof(void*)); - for (i = 0; i < idx->n; ++i) idx->index[i] = kh_init(i); - idx->index2 = (bam_lidx_t*)calloc(idx->n, sizeof(bam_lidx_t)); - - save_bin = save_tid = last_tid = last_bin = 0xffffffffu; - save_off = last_off = bam_tell(fp); last_coor = 0xffffffffu; - n_mapped = n_unmapped = n_no_coor = off_end = 0; - off_beg = off_end = bam_tell(fp); - while ((ret = bam_read1(fp, b)) >= 0) { - if (c->tid < 0) ++n_no_coor; - if (last_tid < c->tid || (last_tid >= 0 && c->tid < 0)) { // change of chromosomes - last_tid = c->tid; - last_bin = 0xffffffffu; - } else if ((uint32_t)last_tid > (uint32_t)c->tid) { - fprintf(stderr, "[bam_index_core] the alignment is not sorted (%s): %d-th chr > %d-th chr\n", - bam1_qname(b), last_tid+1, c->tid+1); - return NULL; - } else if ((int32_t)c->tid >= 0 && last_coor > c->pos) { - fprintf(stderr, "[bam_index_core] the alignment is not sorted (%s): %u > %u in %d-th chr\n", - bam1_qname(b), last_coor, c->pos, c->tid+1); - return NULL; - } - if (c->tid >= 0 && !(c->flag & BAM_FUNMAP)) insert_offset2(&idx->index2[b->core.tid], b, last_off); - if (c->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); - if (last_bin == 0xffffffffu && save_tid != 0xffffffffu) { // write the meta element - off_end = last_off; - insert_offset(idx->index[save_tid], BAM_MAX_BIN, off_beg, off_end); - insert_offset(idx->index[save_tid], BAM_MAX_BIN, n_mapped, n_unmapped); - n_mapped = n_unmapped = 0; - off_beg = off_end; - } - save_off = last_off; - save_bin = last_bin = c->bin; - save_tid = c->tid; - if (save_tid < 0) break; - } - if (bam_tell(fp) <= last_off) { - fprintf(stderr, "[bam_index_core] bug in BGZF/RAZF: %llx < %llx\n", - (unsigned long long)bam_tell(fp), (unsigned long long)last_off); - return NULL; - } - if (c->flag & BAM_FUNMAP) ++n_unmapped; - else ++n_mapped; - last_off = bam_tell(fp); - last_coor = b->core.pos; - } - if (save_tid >= 0) { - insert_offset(idx->index[save_tid], save_bin, save_off, bam_tell(fp)); - insert_offset(idx->index[save_tid], BAM_MAX_BIN, off_beg, bam_tell(fp)); - insert_offset(idx->index[save_tid], BAM_MAX_BIN, n_mapped, n_unmapped); - } - merge_chunks(idx); - fill_missing(idx); - if (ret >= 0) { - while ((ret = bam_read1(fp, b)) >= 0) { - ++n_no_coor; - if (c->tid >= 0 && n_no_coor) { - fprintf(stderr, "[bam_index_core] the alignment is not sorted: reads without coordinates prior to reads with coordinates.\n"); - return NULL; - } - } - } - if (ret < -1) fprintf(stderr, "[bam_index_core] truncated file? Continue anyway. (%d)\n", ret); - free(b->data); free(b); - idx->n_no_coor = n_no_coor; - return idx; -} - -void bam_index_destroy(bam_index_t *idx) -{ - khint_t k; - int i; - if (idx == 0) return; - for (i = 0; i < idx->n; ++i) { - khash_t(i) *index = idx->index[i]; - bam_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); free(idx->index2); - free(idx); -} - -void bam_index_save(const bam_index_t *idx, FILE *fp) -{ - int32_t i, size; - khint_t k; - fwrite("BAI\1", 1, 4, fp); - if (bam_is_be) { - uint32_t x = idx->n; - fwrite(bam_swap_endian_4p(&x), 4, 1, fp); - } else fwrite(&idx->n, 4, 1, fp); - for (i = 0; i < idx->n; ++i) { - khash_t(i) *index = idx->index[i]; - bam_lidx_t *index2 = idx->index2 + i; - // write binning index - size = kh_size(index); - if (bam_is_be) { // big endian - uint32_t x = size; - fwrite(bam_swap_endian_4p(&x), 4, 1, fp); - } else fwrite(&size, 4, 1, fp); - for (k = kh_begin(index); k != kh_end(index); ++k) { - if (kh_exist(index, k)) { - bam_binlist_t *p = &kh_value(index, k); - if (bam_is_be) { // big endian - uint32_t x; - x = kh_key(index, k); fwrite(bam_swap_endian_4p(&x), 4, 1, fp); - x = p->n; fwrite(bam_swap_endian_4p(&x), 4, 1, fp); - for (x = 0; (int)x < p->n; ++x) { - bam_swap_endian_8p(&p->list[x].u); - bam_swap_endian_8p(&p->list[x].v); - } - fwrite(p->list, 16, p->n, fp); - 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 { - fwrite(&kh_key(index, k), 4, 1, fp); - fwrite(&p->n, 4, 1, fp); - fwrite(p->list, 16, p->n, fp); - } - } - } - // write linear index (index2) - if (bam_is_be) { - int x = index2->n; - fwrite(bam_swap_endian_4p(&x), 4, 1, fp); - } else fwrite(&index2->n, 4, 1, fp); - if (bam_is_be) { // big endian - int x; - for (x = 0; (int)x < index2->n; ++x) - bam_swap_endian_8p(&index2->offset[x]); - fwrite(index2->offset, 8, index2->n, fp); - for (x = 0; (int)x < index2->n; ++x) - bam_swap_endian_8p(&index2->offset[x]); - } else fwrite(index2->offset, 8, index2->n, fp); - } - { // write the number of reads coor-less records. - uint64_t x = idx->n_no_coor; - if (bam_is_be) bam_swap_endian_8p(&x); - fwrite(&x, 8, 1, fp); - } - fflush(fp); -} - -static bam_index_t *bam_index_load_core(FILE *fp) -{ - int i; - char magic[4]; - bam_index_t *idx; - if (fp == 0) { - fprintf(stderr, "[bam_index_load_core] fail to load index.\n"); - return 0; - } - fread(magic, 1, 4, fp); - if (strncmp(magic, "BAI\1", 4)) { - fprintf(stderr, "[bam_index_load] wrong magic number.\n"); - fclose(fp); - return 0; - } - idx = (bam_index_t*)calloc(1, sizeof(bam_index_t)); - fread(&idx->n, 4, 1, fp); - if (bam_is_be) bam_swap_endian_4p(&idx->n); - idx->index = (khash_t(i)**)calloc(idx->n, sizeof(void*)); - idx->index2 = (bam_lidx_t*)calloc(idx->n, sizeof(bam_lidx_t)); - for (i = 0; i < idx->n; ++i) { - khash_t(i) *index; - bam_lidx_t *index2 = idx->index2 + i; - uint32_t key, size; - khint_t k; - int j, ret; - bam_binlist_t *p; - index = idx->index[i] = kh_init(i); - // load binning index - fread(&size, 4, 1, fp); - if (bam_is_be) bam_swap_endian_4p(&size); - for (j = 0; j < (int)size; ++j) { - fread(&key, 4, 1, fp); - if (bam_is_be) bam_swap_endian_4p(&key); - k = kh_put(i, index, key, &ret); - p = &kh_value(index, k); - fread(&p->n, 4, 1, fp); - if (bam_is_be) bam_swap_endian_4p(&p->n); - p->m = p->n; - p->list = (pair64_t*)malloc(p->m * 16); - fread(p->list, 16, p->n, fp); - if (bam_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 - fread(&index2->n, 4, 1, fp); - if (bam_is_be) bam_swap_endian_4p(&index2->n); - index2->m = index2->n; - index2->offset = (uint64_t*)calloc(index2->m, 8); - fread(index2->offset, index2->n, 8, fp); - if (bam_is_be) - for (j = 0; j < index2->n; ++j) bam_swap_endian_8p(&index2->offset[j]); - } - if (fread(&idx->n_no_coor, 8, 1, fp) == 0) idx->n_no_coor = 0; - if (bam_is_be) bam_swap_endian_8p(&idx->n_no_coor); - return idx; -} - -bam_index_t *bam_index_load_local(const char *_fn) -{ - FILE *fp; - char *fnidx, *fn; - - if (strstr(_fn, "ftp://") == _fn || strstr(_fn, "http://") == _fn) { - const char *p; - int l = strlen(_fn); - for (p = _fn + l - 1; p >= _fn; --p) - if (*p == '/') break; - fn = strdup(p + 1); - } else fn = strdup(_fn); - fnidx = (char*)calloc(strlen(fn) + 5, 1); - strcpy(fnidx, fn); strcat(fnidx, ".bai"); - fp = fopen(fnidx, "rb"); - if (fp == 0) { // try "{base}.bai" - char *s = strstr(fn, "bam"); - if (s == fn + strlen(fn) - 3) { - strcpy(fnidx, fn); - fnidx[strlen(fn)-1] = 'i'; - fp = fopen(fnidx, "rb"); - } - } - free(fnidx); free(fn); - if (fp) { - bam_index_t *idx = bam_index_load_core(fp); - fclose(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, "wb")) == 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 - -bam_index_t *bam_index_load(const char *fn) -{ - bam_index_t *idx; - idx = bam_index_load_local(fn); - if (idx == 0 && (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn)) { - char *fnidx = calloc(strlen(fn) + 5, 1); - strcat(strcpy(fnidx, fn), ".bai"); - fprintf(stderr, "[bam_index_load] attempting to download the remote index file.\n"); - download_from_remote(fnidx); - free(fnidx); - idx = bam_index_load_local(fn); - } - if (idx == 0) fprintf(stderr, "[bam_index_load] fail to load BAM index.\n"); - return idx; -} - -int bam_index_build2(const char *fn, const char *_fnidx) -{ - char *fnidx; - FILE *fpidx; - bamFile fp; - bam_index_t *idx; - if ((fp = bam_open(fn, "r")) == 0) { - fprintf(stderr, "[bam_index_build2] fail to open the BAM file.\n"); - return -1; - } - idx = bam_index_core(fp); - bam_close(fp); - if(idx == 0) { - fprintf(stderr, "[bam_index_build2] fail to index the BAM file.\n"); - return -1; - } - if (_fnidx == 0) { - fnidx = (char*)calloc(strlen(fn) + 5, 1); - strcpy(fnidx, fn); strcat(fnidx, ".bai"); - } else fnidx = strdup(_fnidx); - fpidx = fopen(fnidx, "wb"); - if (fpidx == 0) { - fprintf(stderr, "[bam_index_build2] fail to create the index file.\n"); - free(fnidx); - bam_index_destroy(idx); - return -1; - } - bam_index_save(idx, fpidx); - bam_index_destroy(idx); - fclose(fpidx); - free(fnidx); - return 0; -} - -int bam_index_build(const char *fn) -{ - return bam_index_build2(fn, 0); -} - -int bam_index(int argc, char *argv[]) -{ - if (argc < 2) { - fprintf(stderr, "Usage: samtools index [out.index]\n"); - return 1; - } - if (argc >= 3) bam_index_build2(argv[1], argv[2]); - else bam_index_build(argv[1]); - return 0; -} - -int bam_idxstats(int argc, char *argv[]) -{ - bam_index_t *idx; - bam_header_t *header; - bamFile fp; - int i; - if (argc < 2) { - fprintf(stderr, "Usage: samtools idxstats \n"); - return 1; - } - fp = bam_open(argv[1], "r"); - if (fp == 0) { fprintf(stderr, "[%s] fail to open BAM.\n", __func__); return 1; } - header = bam_header_read(fp); - bam_close(fp); - idx = bam_index_load(argv[1]); - if (idx == 0) { fprintf(stderr, "[%s] fail to load the index.\n", __func__); return 1; } - for (i = 0; i < idx->n; ++i) { - khint_t k; - khash_t(i) *h = idx->index[i]; - printf("%s\t%d", header->target_name[i], header->target_len[i]); - k = kh_get(i, h, BAM_MAX_BIN); - if (k != kh_end(h)) - printf("\t%llu\t%llu", (long long)kh_val(h, k).list[1].u, (long long)kh_val(h, k).list[1].v); - else printf("\t0\t0"); - putchar('\n'); - } - printf("*\t0\t0\t%llu\n", (long long)idx->n_no_coor); - bam_header_destroy(header); - bam_index_destroy(idx); - return 0; -} - -static inline int reg2bins(uint32_t beg, uint32_t end, uint16_t list[BAM_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; -} - -static inline int is_overlap(uint32_t beg, uint32_t end, const bam1_t *b) -{ - uint32_t rbeg = b->core.pos; - uint32_t rend = b->core.n_cigar? bam_calend(&b->core, bam1_cigar(b)) : b->core.pos + 1; - return (rend > beg && rbeg < end); -} - -struct __bam_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; - pair64_t *off; -}; - -// bam_fetch helper function retrieves -bam_iter_t bam_iter_query(const bam_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; - bam_iter_t iter = 0; - - if (beg < 0) beg = 0; - if (end < beg) return 0; - // initialize iter - iter = calloc(1, sizeof(struct __bam_iter_t)); - iter->tid = tid, iter->beg = beg, iter->end = end; iter->i = -1; - // - bins = (uint16_t*)calloc(BAM_MAX_BIN, 2); - n_bins = reg2bins(beg, end, bins); - index = idx->index[tid]; - if (idx->index2[tid].n > 0) { - min_off = (beg>>BAM_LIDX_SHIFT >= idx->index2[tid].n)? idx->index2[tid].offset[idx->index2[tid].n-1] - : idx->index2[tid].offset[beg>>BAM_LIDX_SHIFT]; - if (min_off == 0) { // improvement for index files built by tabix prior to 0.1.4 - int n = beg>>BAM_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; - bam_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]; - } - } - free(bins); - if (n_off == 0) { - free(off); return iter; - } - { - bam1_t *b = (bam1_t*)calloc(1, sizeof(bam1_t)); - int l; - ks_introsort(off, 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 -#if defined(BAM_TRUE_OFFSET) || defined(BAM_VIRTUAL_OFFSET16) - for (i = 1, l = 0; i < n_off; ++i) { -#ifdef BAM_TRUE_OFFSET - if (off[l].v + BAM_MIN_CHUNK_GAP > off[i].u) off[l].v = off[i].v; -#else - if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v; -#endif - else off[++l] = off[i]; - } - n_off = l + 1; -#endif - } - bam_destroy1(b); - } - iter->n_off = n_off; iter->off = off; - return iter; -} - -pair64_t *get_chunk_coordinates(const bam_index_t *idx, int tid, int beg, int end, int *cnt_off) -{ // for pysam compatibility - bam_iter_t iter; - pair64_t *off; - iter = bam_iter_query(idx, tid, beg, end); - off = iter->off; *cnt_off = iter->n_off; - free(iter); - return off; -} - -void bam_iter_destroy(bam_iter_t iter) -{ - if (iter) { free(iter->off); free(iter); } -} - -int bam_iter_read(bamFile fp, bam_iter_t iter, bam1_t *b) -{ - int ret; - if (iter && iter->finished) return -1; - if (iter == 0 || iter->from_first) { - ret = bam_read1(fp, b); - if (ret < 0 && iter) iter->finished = 1; - return ret; - } - if (iter->off == 0) return -1; - for (;;) { - 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) { ret = -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 - bam_seek(fp, iter->off[iter->i+1].u, SEEK_SET); - iter->curr_off = bam_tell(fp); - } - ++iter->i; - } - if ((ret = bam_read1(fp, b)) >= 0) { - iter->curr_off = bam_tell(fp); - if (b->core.tid != iter->tid || b->core.pos >= iter->end) { // no need to proceed - ret = bam_validate1(NULL, b)? -1 : -5; // determine whether end of region or error - break; - } - else if (is_overlap(iter->beg, iter->end, b)) return ret; - } else break; // end of file or error - } - iter->finished = 1; - return ret; -} - -int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func) -{ - int ret; - bam_iter_t iter; - bam1_t *b; - b = bam_init1(); - iter = bam_iter_query(idx, tid, beg, end); - while ((ret = bam_iter_read(fp, iter, b)) >= 0) func(b, data); - bam_iter_destroy(iter); - bam_destroy1(b); - return (ret == -1)? 0 : ret; -} diff --git a/samlib/sam_header.c b/samlib/sam_header.c deleted file mode 100644 index 88b6a1c..0000000 --- a/samlib/sam_header.c +++ /dev/null @@ -1,810 +0,0 @@ -#include "sam_header.h" -#include -#include -#include -#include -#include - -#include "khash.h" -KHASH_MAP_INIT_STR(str, const char *) - -struct _HeaderList -{ - struct _HeaderList *last; // Hack: Used and maintained only by list_append_to_end. Maintained in the root node only. - struct _HeaderList *next; - void *data; -}; -typedef struct _HeaderList list_t; -typedef list_t HeaderDict; - -typedef struct -{ - char key[2]; - char *value; -} -HeaderTag; - -typedef struct -{ - char type[2]; - list_t *tags; -} -HeaderLine; - -const char *o_hd_tags[] = {"SO","GO",NULL}; -const char *r_hd_tags[] = {"VN",NULL}; - -const char *o_sq_tags[] = {"AS","M5","UR","SP",NULL}; -const char *r_sq_tags[] = {"SN","LN",NULL}; -const char *u_sq_tags[] = {"SN",NULL}; - -const char *o_rg_tags[] = {"CN","DS","DT","FO","KS","LB","PG","PI","PL","PU","SM",NULL}; -const char *r_rg_tags[] = {"ID",NULL}; -const char *u_rg_tags[] = {"ID",NULL}; - -const char *o_pg_tags[] = {"VN","CL",NULL}; -const char *r_pg_tags[] = {"ID",NULL}; - -const char *types[] = {"HD","SQ","RG","PG","CO",NULL}; -const char **optional_tags[] = {o_hd_tags,o_sq_tags,o_rg_tags,o_pg_tags,NULL,NULL}; -const char **required_tags[] = {r_hd_tags,r_sq_tags,r_rg_tags,r_pg_tags,NULL,NULL}; -const char **unique_tags[] = {NULL, u_sq_tags,u_rg_tags,NULL,NULL,NULL}; - - -static void debug(const char *format, ...) -{ - va_list ap; - va_start(ap, format); - vfprintf(stderr, format, ap); - va_end(ap); -} - -#if 0 -// Replaced by list_append_to_end -static list_t *list_prepend(list_t *root, void *data) -{ - list_t *l = malloc(sizeof(list_t)); - l->next = root; - l->data = data; - return l; -} -#endif - -// Relies on the root->last being correct. Do not use with the other list_* -// routines unless they are fixed to modify root->last as well. -static list_t *list_append_to_end(list_t *root, void *data) -{ - list_t *l = malloc(sizeof(list_t)); - l->last = l; - l->next = NULL; - l->data = data; - - if ( !root ) - return l; - - root->last->next = l; - root->last = l; - return root; -} - -static list_t *list_append(list_t *root, void *data) -{ - list_t *l = root; - while (l && l->next) - l = l->next; - if ( l ) - { - l->next = malloc(sizeof(list_t)); - l = l->next; - } - else - { - l = malloc(sizeof(list_t)); - root = l; - } - l->data = data; - l->next = NULL; - return root; -} - -static void list_free(list_t *root) -{ - list_t *l = root; - while (root) - { - l = root; - root = root->next; - free(l); - } -} - - - -// Look for a tag "XY" in a predefined const char *[] array. -static int tag_exists(const char *tag, const char **tags) -{ - int itag=0; - if ( !tags ) return -1; - while ( tags[itag] ) - { - if ( tags[itag][0]==tag[0] && tags[itag][1]==tag[1] ) return itag; - itag++; - } - return -1; -} - - - -// Mimics the behaviour of getline, except it returns pointer to the next chunk of the text -// or NULL if everything has been read. The lineptr should be freed by the caller. The -// newline character is stripped. -static const char *nextline(char **lineptr, size_t *n, const char *text) -{ - int len; - const char *to = text; - - if ( !*to ) return NULL; - - while ( *to && *to!='\n' && *to!='\r' ) to++; - len = to - text + 1; - - if ( *to ) - { - // Advance the pointer for the next call - if ( *to=='\n' ) to++; - else if ( *to=='\r' && *(to+1)=='\n' ) to+=2; - } - if ( !len ) - return to; - - if ( !*lineptr ) - { - *lineptr = malloc(len); - *n = len; - } - else if ( *nkey[0] = name[0]; - tag->key[1] = name[1]; - tag->value = malloc(len+1); - memcpy(tag->value,value_from,len+1); - tag->value[len] = 0; - return tag; -} - -static HeaderTag *header_line_has_tag(HeaderLine *hline, const char *key) -{ - list_t *tags = hline->tags; - while (tags) - { - HeaderTag *tag = tags->data; - if ( tag->key[0]==key[0] && tag->key[1]==key[1] ) return tag; - tags = tags->next; - } - return NULL; -} - - -// Return codes: -// 0 .. different types or unique tags differ or conflicting tags, cannot be merged -// 1 .. all tags identical -> no need to merge, drop one -// 2 .. the unique tags match and there are some conflicting tags (same tag, different value) -> error, cannot be merged nor duplicated -// 3 .. there are some missing complementary tags and no unique conflict -> can be merged into a single line -static int sam_header_compare_lines(HeaderLine *hline1, HeaderLine *hline2) -{ - HeaderTag *t1, *t2; - - if ( hline1->type[0]!=hline2->type[0] || hline1->type[1]!=hline2->type[1] ) - return 0; - - int itype = tag_exists(hline1->type,types); - if ( itype==-1 ) { - debug("[sam_header_compare_lines] Unknown type [%c%c]\n", hline1->type[0],hline1->type[1]); - return -1; // FIXME (lh3): error; I do not know how this will be handled in Petr's code - } - - if ( unique_tags[itype] ) - { - t1 = header_line_has_tag(hline1,unique_tags[itype][0]); - t2 = header_line_has_tag(hline2,unique_tags[itype][0]); - if ( !t1 || !t2 ) // this should never happen, the unique tags are required - return 2; - - if ( strcmp(t1->value,t2->value) ) - return 0; // the unique tags differ, cannot be merged - } - if ( !required_tags[itype] && !optional_tags[itype] ) - { - t1 = hline1->tags->data; - t2 = hline2->tags->data; - if ( !strcmp(t1->value,t2->value) ) return 1; // identical comments - return 0; - } - - int missing=0, itag=0; - while ( required_tags[itype] && required_tags[itype][itag] ) - { - t1 = header_line_has_tag(hline1,required_tags[itype][itag]); - t2 = header_line_has_tag(hline2,required_tags[itype][itag]); - if ( !t1 && !t2 ) - return 2; // this should never happen - else if ( !t1 || !t2 ) - missing = 1; // there is some tag missing in one of the hlines - else if ( strcmp(t1->value,t2->value) ) - { - if ( unique_tags[itype] ) - return 2; // the lines have a matching unique tag but have a conflicting tag - - return 0; // the lines contain conflicting tags, cannot be merged - } - itag++; - } - itag = 0; - while ( optional_tags[itype] && optional_tags[itype][itag] ) - { - t1 = header_line_has_tag(hline1,optional_tags[itype][itag]); - t2 = header_line_has_tag(hline2,optional_tags[itype][itag]); - if ( !t1 && !t2 ) - { - itag++; - continue; - } - if ( !t1 || !t2 ) - missing = 1; // there is some tag missing in one of the hlines - else if ( strcmp(t1->value,t2->value) ) - { - if ( unique_tags[itype] ) - return 2; // the lines have a matching unique tag but have a conflicting tag - - return 0; // the lines contain conflicting tags, cannot be merged - } - itag++; - } - if ( missing ) return 3; // there are some missing complementary tags with no conflicts, can be merged - return 1; -} - - -static HeaderLine *sam_header_line_clone(const HeaderLine *hline) -{ - list_t *tags; - HeaderLine *out = malloc(sizeof(HeaderLine)); - out->type[0] = hline->type[0]; - out->type[1] = hline->type[1]; - out->tags = NULL; - - tags = hline->tags; - while (tags) - { - HeaderTag *old = tags->data; - - HeaderTag *new = malloc(sizeof(HeaderTag)); - new->key[0] = old->key[0]; - new->key[1] = old->key[1]; - new->value = strdup(old->value); - out->tags = list_append(out->tags, new); - - tags = tags->next; - } - return out; -} - -static int sam_header_line_merge_with(HeaderLine *out_hline, const HeaderLine *tmpl_hline) -{ - list_t *tmpl_tags; - - if ( out_hline->type[0]!=tmpl_hline->type[0] || out_hline->type[1]!=tmpl_hline->type[1] ) - return 0; - - tmpl_tags = tmpl_hline->tags; - while (tmpl_tags) - { - HeaderTag *tmpl_tag = tmpl_tags->data; - HeaderTag *out_tag = header_line_has_tag(out_hline, tmpl_tag->key); - if ( !out_tag ) - { - HeaderTag *tag = malloc(sizeof(HeaderTag)); - tag->key[0] = tmpl_tag->key[0]; - tag->key[1] = tmpl_tag->key[1]; - tag->value = strdup(tmpl_tag->value); - out_hline->tags = list_append(out_hline->tags,tag); - } - tmpl_tags = tmpl_tags->next; - } - return 1; -} - - -static HeaderLine *sam_header_line_parse(const char *headerLine) -{ - HeaderLine *hline; - HeaderTag *tag; - const char *from, *to; - from = headerLine; - - if ( *from != '@' ) { - debug("[sam_header_line_parse] expected '@', got [%s]\n", headerLine); - return 0; - } - to = ++from; - - while (*to && *to!='\t') to++; - if ( to-from != 2 ) { - debug("[sam_header_line_parse] expected '@XY', got [%s]\nHint: The header tags must be tab-separated.\n", headerLine); - return 0; - } - - hline = malloc(sizeof(HeaderLine)); - hline->type[0] = from[0]; - hline->type[1] = from[1]; - hline->tags = NULL; - - int itype = tag_exists(hline->type, types); - - from = to; - while (*to && *to=='\t') to++; - if ( to-from != 1 ) { - debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from)); - free(hline); - return 0; - } - from = to; - while (*from) - { - while (*to && *to!='\t') to++; - - if ( !required_tags[itype] && !optional_tags[itype] ) - { - // CO is a special case, it can contain anything, including tabs - if ( *to ) { to++; continue; } - tag = new_tag(" ",from,to-1); - } - else - tag = new_tag(from,from+3,to-1); - - if ( header_line_has_tag(hline,tag->key) ) - debug("The tag '%c%c' present (at least) twice on line [%s]\n", tag->key[0],tag->key[1], headerLine); - hline->tags = list_append(hline->tags, tag); - - from = to; - while (*to && *to=='\t') to++; - if ( *to && to-from != 1 ) { - debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from)); - return 0; - } - - from = to; - } - return hline; -} - - -// Must be of an existing type, all tags must be recognised and all required tags must be present -static int sam_header_line_validate(HeaderLine *hline) -{ - list_t *tags; - HeaderTag *tag; - int itype, itag; - - // Is the type correct? - itype = tag_exists(hline->type, types); - if ( itype==-1 ) - { - debug("The type [%c%c] not recognised.\n", hline->type[0],hline->type[1]); - return 0; - } - - // Has all required tags? - itag = 0; - while ( required_tags[itype] && required_tags[itype][itag] ) - { - if ( !header_line_has_tag(hline,required_tags[itype][itag]) ) - { - debug("The tag [%c%c] required for [%c%c] not present.\n", required_tags[itype][itag][0],required_tags[itype][itag][1], - hline->type[0],hline->type[1]); - return 0; - } - itag++; - } - - // Are all tags recognised? - tags = hline->tags; - while ( tags ) - { - tag = tags->data; - if ( !tag_exists(tag->key,required_tags[itype]) && !tag_exists(tag->key,optional_tags[itype]) ) - { - // Lower case tags are user-defined values. - if( !(islower(tag->key[0]) || islower(tag->key[1])) ) - { - // Neither is lower case, but tag was not recognized. - debug("Unknown tag [%c%c] for [%c%c].\n", tag->key[0],tag->key[1], hline->type[0],hline->type[1]); - // return 0; // Even unknown tags are allowed - for forward compatibility with new attributes - } - // else - allow user defined tag - } - tags = tags->next; - } - - return 1; -} - - -static void print_header_line(FILE *fp, HeaderLine *hline) -{ - list_t *tags = hline->tags; - HeaderTag *tag; - - fprintf(fp, "@%c%c", hline->type[0],hline->type[1]); - while (tags) - { - tag = tags->data; - - fprintf(fp, "\t"); - if ( tag->key[0]!=' ' || tag->key[1]!=' ' ) - fprintf(fp, "%c%c:", tag->key[0],tag->key[1]); - fprintf(fp, "%s", tag->value); - - tags = tags->next; - } - fprintf(fp,"\n"); -} - - -static void sam_header_line_free(HeaderLine *hline) -{ - list_t *tags = hline->tags; - while (tags) - { - HeaderTag *tag = tags->data; - free(tag->value); - free(tag); - tags = tags->next; - } - list_free(hline->tags); - free(hline); -} - -void sam_header_free(void *_header) -{ - HeaderDict *header = (HeaderDict*)_header; - list_t *hlines = header; - while (hlines) - { - sam_header_line_free(hlines->data); - hlines = hlines->next; - } - list_free(header); -} - -HeaderDict *sam_header_clone(const HeaderDict *dict) -{ - HeaderDict *out = NULL; - while (dict) - { - HeaderLine *hline = dict->data; - out = list_append(out, sam_header_line_clone(hline)); - dict = dict->next; - } - return out; -} - -// Returns a newly allocated string -char *sam_header_write(const void *_header) -{ - const HeaderDict *header = (const HeaderDict*)_header; - char *out = NULL; - int len=0, nout=0; - const list_t *hlines; - - // Calculate the length of the string to allocate - hlines = header; - while (hlines) - { - len += 4; // @XY and \n - - HeaderLine *hline = hlines->data; - list_t *tags = hline->tags; - while (tags) - { - HeaderTag *tag = tags->data; - len += strlen(tag->value) + 1; // \t - if ( tag->key[0]!=' ' || tag->key[1]!=' ' ) - len += strlen(tag->value) + 3; // XY: - tags = tags->next; - } - hlines = hlines->next; - } - - nout = 0; - out = malloc(len+1); - hlines = header; - while (hlines) - { - HeaderLine *hline = hlines->data; - - nout += sprintf(out+nout,"@%c%c",hline->type[0],hline->type[1]); - - list_t *tags = hline->tags; - while (tags) - { - HeaderTag *tag = tags->data; - nout += sprintf(out+nout,"\t"); - if ( tag->key[0]!=' ' || tag->key[1]!=' ' ) - nout += sprintf(out+nout,"%c%c:", tag->key[0],tag->key[1]); - nout += sprintf(out+nout,"%s", tag->value); - tags = tags->next; - } - hlines = hlines->next; - nout += sprintf(out+nout,"\n"); - } - out[len] = 0; - return out; -} - -void *sam_header_parse2(const char *headerText) -{ - list_t *hlines = NULL; - HeaderLine *hline; - const char *text; - char *buf=NULL; - size_t nbuf = 0; - int tovalidate = 0; - - if ( !headerText ) - return 0; - - text = headerText; - while ( (text=nextline(&buf, &nbuf, text)) ) - { - hline = sam_header_line_parse(buf); - if ( hline && (!tovalidate || sam_header_line_validate(hline)) ) - // With too many (~250,000) reference sequences the header parsing was too slow with list_append. - hlines = list_append_to_end(hlines, hline); - else - { - if (hline) sam_header_line_free(hline); - sam_header_free(hlines); - if ( buf ) free(buf); - return NULL; - } - } - if ( buf ) free(buf); - - return hlines; -} - -void *sam_header2tbl(const void *_dict, char type[2], char key_tag[2], char value_tag[2]) -{ - const HeaderDict *dict = (const HeaderDict*)_dict; - const list_t *l = dict; - khash_t(str) *tbl = kh_init(str); - khiter_t k; - int ret; - - if (_dict == 0) return tbl; // return an empty (not null) hash table - while (l) - { - HeaderLine *hline = l->data; - if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] ) - { - l = l->next; - continue; - } - - HeaderTag *key, *value; - key = header_line_has_tag(hline,key_tag); - value = header_line_has_tag(hline,value_tag); - if ( !key || !value ) - { - l = l->next; - continue; - } - - k = kh_get(str, tbl, key->value); - if ( k != kh_end(tbl) ) - debug("[sam_header_lookup_table] They key %s not unique.\n", key->value); - k = kh_put(str, tbl, key->value, &ret); - kh_value(tbl, k) = value->value; - - l = l->next; - } - return tbl; -} - -char **sam_header2list(const void *_dict, char type[2], char key_tag[2], int *_n) -{ - const HeaderDict *dict = (const HeaderDict*)_dict; - const list_t *l = dict; - int max, n; - char **ret; - - ret = 0; *_n = max = n = 0; - while (l) - { - HeaderLine *hline = l->data; - if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] ) - { - l = l->next; - continue; - } - - HeaderTag *key; - key = header_line_has_tag(hline,key_tag); - if ( !key ) - { - l = l->next; - continue; - } - - if (n == max) { - max = max? max<<1 : 4; - ret = realloc(ret, max * sizeof(void*)); - } - ret[n++] = key->value; - - l = l->next; - } - *_n = n; - return ret; -} - -void *sam_header2key_val(void *iter, const char type[2], const char key_tag[2], const char value_tag[2], const char **_key, const char **_value) -{ - list_t *l = iter; - if ( !l ) return NULL; - - while (l) - { - HeaderLine *hline = l->data; - if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] ) - { - l = l->next; - continue; - } - - HeaderTag *key, *value; - key = header_line_has_tag(hline,key_tag); - value = header_line_has_tag(hline,value_tag); - if ( !key && !value ) - { - l = l->next; - continue; - } - - *_key = key->value; - *_value = value->value; - return l->next; - } - return l; -} - -const char *sam_tbl_get(void *h, const char *key) -{ - khash_t(str) *tbl = (khash_t(str)*)h; - khint_t k; - k = kh_get(str, tbl, key); - return k == kh_end(tbl)? 0 : kh_val(tbl, k); -} - -int sam_tbl_size(void *h) -{ - khash_t(str) *tbl = (khash_t(str)*)h; - return h? kh_size(tbl) : 0; -} - -void sam_tbl_destroy(void *h) -{ - khash_t(str) *tbl = (khash_t(str)*)h; - kh_destroy(str, tbl); -} - -void *sam_header_merge(int n, const void **_dicts) -{ - const HeaderDict **dicts = (const HeaderDict**)_dicts; - HeaderDict *out_dict; - int idict, status; - - if ( n<2 ) return NULL; - - out_dict = sam_header_clone(dicts[0]); - - for (idict=1; idictdata, out_hlines->data); - if ( status==0 ) - { - out_hlines = out_hlines->next; - continue; - } - - if ( status==2 ) - { - print_header_line(stderr,tmpl_hlines->data); - print_header_line(stderr,out_hlines->data); - debug("Conflicting lines, cannot merge the headers.\n"); - return 0; - } - if ( status==3 ) - sam_header_line_merge_with(out_hlines->data, tmpl_hlines->data); - - inserted = 1; - break; - } - if ( !inserted ) - out_dict = list_append(out_dict, sam_header_line_clone(tmpl_hlines->data)); - - tmpl_hlines = tmpl_hlines->next; - } - } - - return out_dict; -} - -char **sam_header2tbl_n(const void *dict, const char type[2], const char *tags[], int *n) -{ - int nout = 0; - char **out = NULL; - - *n = 0; - list_t *l = (list_t *)dict; - if ( !l ) return NULL; - - int i, ntags = 0; - while ( tags[ntags] ) ntags++; - - while (l) - { - HeaderLine *hline = l->data; - if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] ) - { - l = l->next; - continue; - } - out = (char**) realloc(out, sizeof(char*)*(nout+1)*ntags); - for (i=0; ivalue; - } - nout++; - l = l->next; - } - *n = nout; - return out; -} - diff --git a/samlib/sam_header.h b/samlib/sam_header.h deleted file mode 100644 index 7f88a21..0000000 --- a/samlib/sam_header.h +++ /dev/null @@ -1,48 +0,0 @@ -#ifndef __SAM_HEADER_H__ -#define __SAM_HEADER_H__ - -#ifdef __cplusplus -extern "C" { -#endif - - void *sam_header_parse2(const char *headerText); - void *sam_header_merge(int n, const void **dicts); - void sam_header_free(void *header); - char *sam_header_write(const void *headerDict); // returns a newly allocated string - - /* - // Usage example - const char *key, *val; - void *iter = sam_header_parse2(bam->header->text); - while ( iter = sam_header_key_val(iter, "RG","ID","SM" &key,&val) ) printf("%s\t%s\n", key,val); - */ - void *sam_header2key_val(void *iter, const char type[2], const char key_tag[2], const char value_tag[2], const char **key, const char **value); - char **sam_header2list(const void *_dict, char type[2], char key_tag[2], int *_n); - - /* - // Usage example - int i, j, n; - const char *tags[] = {"SN","LN","UR","M5",NULL}; - void *dict = sam_header_parse2(bam->header->text); - char **tbl = sam_header2tbl_n(h->dict, "SQ", tags, &n); - for (i=0; i + +// htslib for BAM/CRAM/SAM file reading +#include +#include +#include // khash, kstring and knetfile are standard utils of klib -#include "bgzf.h" // write tabix-able depth.gz file #include "khash.h" #include "knetfile.h" #include "kstring.h" #include -static char const *program_name = "bamdst"; -static char const *Version = "1.1.0"; +static char const *program_name = "xamdst"; +static char const *Version = "2.0.0"; /* flank region will be stat in the coverage report file, * this value can be set by -f / --flank */ @@ -57,9 +61,11 @@ static bool stdin_lock = FALSE; static bool zero_based = TRUE; /* 定义 max cutoff 最大值为10 */ - static const int MAX_CUTOFFS = 10; +/* I/O 缓冲区大小优化 - 64KB 缓冲区减少系统调用次数 */ +static const int WRITE_BUFFER_SIZE = 65536; + /* The number of threads after which there are diminishing performance gains. */ // enum { DEFAULT_MAX_THREADS = 8 }; @@ -71,7 +77,7 @@ static const int MAX_CUTOFFS = 10; /* export target reads to a specitified bam file */ static char *export_target_bam = NULL; -bamFile bamoutfp; +htsFile *bamoutfp = NULL; int check_filename_isbam(char *name) { @@ -103,21 +109,19 @@ void h_chrlength_init() h_chrlen = kh_init(chr); } -// warp bamheader to retrieve chromosome legth hash -void header2chrhash(bam_header_t *h) +// warp sam header to retrieve chromosome length hash +void header2chrhash(sam_hdr_t *h) { - int i, n, ret; + int i, ret; khiter_t k; - h->dict = sam_header_parse2(h->text); - const char *tags[] = {"SN", "LN", "UR", "M5", NULL}; - char **tbl = sam_header2tbl_n(h->dict, "SQ", tags, &n); - for (i = 0; i < n; i++) + int n_targets = sam_hdr_nref(h); + for (i = 0; i < n_targets; i++) { - k = kh_put(chr, h_chrlen, strdup(tbl[4 * i]), &ret); - kh_val(h_chrlen, k).length = atoi(tbl[4 * i + 1]); + const char *name = sam_hdr_tid2name(h, i); + hts_pos_t len = sam_hdr_tid2len(h, i); + k = kh_put(chr, h_chrlen, strdup(name), &ret); + kh_val(h_chrlen, k).length = (uint32_t)len; } - if (tbl) - free(tbl); } void chrhash_destroy() @@ -154,6 +158,15 @@ struct opt_aux int *cutoffs; // 指向 cutoff 数组的指针 int mapQ_lim; int maxdepth; + // 新增：深度比例参数 + bool depth_ratio; // 是否启用深度比例统计 + int num_ratios; // 比例数量 + int max_ratios; // 最大比例数量 + float *ratios; // 比例数组 (如 0.1, 0.2, 0.5) + // 新增：参考基因组路径（用于 CRAM 文件） + char *reference; // 参考基因组 FASTA 文件路径 + // 新增：线程数 + int nthreads; // htslib 多线程压缩/解压缩线程数 }; // 一个函数来初始化 opt_aux 结构体 @@ -175,6 +188,27 @@ struct opt_aux init_opt_aux() opt.inputs = NULL; opt.isize_lim = 2000; opt.mapQ_lim = 20; + + // 初始化深度比例参数 + opt.depth_ratio = TRUE; // 默认启用 + opt.num_ratios = 2; // 默认2个比例 + opt.max_ratios = MAX_CUTOFFS; + opt.ratios = malloc(opt.max_ratios * sizeof(float)); + if (opt.ratios == NULL) + { + fprintf(stderr, "Memory allocation failed for ratios\n"); + exit(EXIT_FAILURE); + } + // 设置默认比例值 0.2 和 0.5 + opt.ratios[0] = 0.2f; + opt.ratios[1] = 0.5f; + + // 初始化参考基因组路径 + opt.reference = NULL; + + // 初始化线程数（0 = 单线程模式，不使用多线程） + opt.nthreads = 0; + return opt; } @@ -184,7 +218,6 @@ struct depnode unsigned start; unsigned stop; unsigned *vals; // raw depth - unsigned *cnts; // reads count in this position unsigned *rmdupdep; // clean depth, rmdup, mapQ > 20, primary hit unsigned *covdep; // coverage depth struct depnode *next; @@ -206,7 +239,6 @@ static const char *init_debugmsg[] = {"Success", "Trying to allocated an unempty } while (0) /* node->vals is the depth value of each loc - * node->cnts is the count of covered reads * init the memory before use it */ static int depnode_init(struct depnode *node) { @@ -218,11 +250,9 @@ static int depnode_init(struct depnode *node) if (isZero(node->len)) return 2; node->vals = (unsigned *)needmem((node->len) * sizeof(unsigned)); - node->cnts = (unsigned *)needmem((node->len) * sizeof(unsigned)); node->rmdupdep = (unsigned *)needmem((node->len) * sizeof(unsigned)); node->covdep = (unsigned *)needmem((node->len) * sizeof(unsigned)); memset(node->vals, 0, node->len * sizeof(unsigned)); - memset(node->cnts, 0, node->len * sizeof(unsigned)); memset(node->rmdupdep, 0, node->len * sizeof(unsigned)); memset(node->covdep, 0, node->len * sizeof(unsigned)); return 0; @@ -237,7 +267,6 @@ static int depnode_init(struct depnode *node) struct depnode *tmpnode = node; \ node = node->next; \ freemem(tmpnode->vals); \ - freemem(tmpnode->cnts); \ freemem(tmpnode->rmdupdep); \ freemem(tmpnode->covdep); \ freemem(tmpnode); \ @@ -261,7 +290,7 @@ static struct depnode *bed_depnode_list(bedreglist_t *bed) // same beg and end pos /* the length of this region should be zero if not allocated memory yet - * Assign the length value when init the vals and cnts */ + * Assign the length value when init the vals */ node->len = 0; if (isZero(i)) header = node; @@ -287,10 +316,10 @@ struct _aux uint64_t tgt_len, flk_len; unsigned tgt_nreg; - /* data, array of bam struct - * h, point to bam header */ - bamFile *data; - bam_header_t *h; + /* data, array of htsFile pointers for BAM/CRAM/SAM + * h, point to sam header */ + htsFile **data; + sam_hdr_t *h; /* h_tgt, target bed hash * h_flk, flank bed hash */ @@ -312,8 +341,8 @@ struct _aux *aux_init() struct _aux *a; a = calloc(1, sizeof(struct _aux)); a->nchr = a->maxdep = a->ndata = 0; - a->data = NULL; // bamFile - a->h = NULL; // bam header + a->data = NULL; // htsFile array + a->h = NULL; // sam header a->h_tgt = kh_init(reg); a->h_flk = kh_init(reg); count32_init(a->c_dep); @@ -335,7 +364,7 @@ void aux_destroy(struct _aux *a) free(a->data); bedHand->destroy((void *)a->h_tgt, destroy_data); bedHand->destroy((void *)a->h_flk, destroy_void); - bam_header_destroy(a->h); + sam_hdr_destroy(a->h); /* if (a->c_dep->n > 0) count_destroy(a->c_dep); */ /* if (a->c_rmdupdep->n > 0) count_destroy(a->c_rmdupdep); */ /* if (a->c_flkdep->n > 0) count_destroy(a->c_flkdep); */ @@ -443,34 +472,43 @@ void usage(int status) printf("\n\ bamdst version: %s\n\ USAGE : %s [OPTION] -p -o [in1.bam [in2.bam ... ]]\n\ + or : %s [OPTION] -p -o [in1.cram] -T \n\ or : %s [OPTION] -p -o -\n\ ", - Version, program_name, program_name); + Version, program_name, program_name, program_name); puts("\ Option -o and -p are mandatory:\n\ - -o, --outdir output dir\n\ + -o, --outdir output dir (will be created if not exists)\n\ -p, --bed probe or target regions file, the region file will \n\ be merged before calculate depths\n\ "); puts("\ Optional parameters:\n\ + -T, --reference FILE reference genome FASTA file (required for CRAM input)\n\ -f, --flank [200] flank n bp of each region\n\ -q [20] map quality cutoff value, greater or equal to the value will be count\n\ --maxdepth [0] set the max depth to stat the cumu distribution.\n\ --cutoffdepth [0,0] list the coverage of above these depths, allow maximal 10 cutoffs.\n\ + --depthratio [0.2,0.5] coverage at ratios of average depth (e.g., 0.1,0.2,0.5)\n\ --isize [2000] stat the inferred insert size under this value\n\ --uncover [5] region will included in uncover file if below it\n\ --bamout BAMFILE target reads will be exported to this bam file\n\ + --threads [0] number of threads for BAM/CRAM I/O (0 = single-threaded)\n\ -1 begin position of bed file is 1-based\n\ -h, --help print this help info\n\ \n"); puts("\ +* Supported input formats: BAM, CRAM, SAM\n\ +* CRAM files require a reference genome (-T/--reference option)\n\ +* The output directory will be created automatically if it doesn't exist\n\ +\n\ * Five essential files would be created in the output dir. \n\ * region.tsv.gz and depth.tsv.gz are zipped by bgzip, so you can use tabix \n\ index these files.\n\n\ - coverage.report a report of the coverage information and reads \n\ information of whole target regions\n\ + - coverage.report.json JSON format of coverage report for easy parsing\n\ - cumu.plot distribution data of depth values\n\ - insert.plot distribution data of inferred insert size \n\ - chromosome.report coverage information for each chromosome\n\ @@ -491,7 +529,7 @@ Optional parameters:\n\ "); puts("============\n"); puts(" HOMEPAGE: \n\ - https://github.com/shiquan/bamdst\n"); + https://github.com/pzweuj/xamdst\n"); } exit(EXIT_SUCCESS); } @@ -549,7 +587,9 @@ static float median_cal(const uint32_t *array, int l) tmp = (uint32_t *)needmem(l * sizeof(uint32_t)); memcpy(tmp, array, l * sizeof(uint32_t)); ks_introsort(uint32_t, l, tmp); - float med = l & 1 ? tmp[(l >> 1) + 1] : (float)(tmp[l >> 1] + tmp[(l >> 1) - 1]) / 2; + // For odd l the median is the middle element tmp[l>>1]; for even l it is + // the average of the two middle elements tmp[l>>1] and tmp[(l>>1)-1]. + float med = l & 1 ? tmp[l >> 1] : (float)(tmp[l >> 1] + tmp[(l >> 1) - 1]) / 2; mustfree(tmp); return med; } @@ -579,11 +619,18 @@ static float coverage_cal(const uint32_t *array, int l) return cov * 100; } -// FIXME: need broken when bed file is truncated +// 加载并初始化 BED 文件，包含截断检测 int load_bed_init(char const *fn, aux_t *a) { int ret = 0; bedHand->read(fn, a->h_tgt, 0, 0, &ret); + + // 检查 BED 文件是否为空 + if (kh_size(a->h_tgt) == 0) + { + errabort("Failed to read BED file: %s (file may be empty or corrupted)", fn); + } + if (zero_based && ret) { warnings("This region is not a standard bed format.\n" @@ -646,17 +693,18 @@ typedef enum UNKNOWN } cntstat_t; -int match_pos(struct depnode *header, uint32_t pos, cntstat_t state) +int match_pos(struct depnode **pp, uint32_t pos, cntstat_t state) { - struct depnode *tmp = header; + struct depnode *tmp = *pp; + /* forward-only advancement: pos grows monotonically along a read, so the + * node pointer only ever moves forward, avoiding a re-walk from the head + * on every base. */ while (tmp && pos > tmp->stop) - { tmp = tmp->next; - } + *pp = tmp; if (isNull(tmp)) return 1; // this chromosome is finished, skip in next loop - // debug("pos: %u\tstart: %u\tstop: %u", pos, tmp->start, tmp->stop); if (pos >= tmp->start) { if (isZero(tmp->len)) @@ -678,7 +726,7 @@ int match_pos(struct depnode *header, uint32_t pos, cntstat_t state) } return 0; } - return 1; // not reachable + return 1; // pos falls in a gap before this region } /* when deal with a read struct, check the begin of this read and the last @@ -699,8 +747,8 @@ int readcore(struct depnode *header, bam1_t const *b, cntstat_t state) /* header = tmp; */ if (isNull(tmp)) return 0; - uint32_t *cigar = bam1_cigar(b); - uint32_t end = bam_calend(c, cigar); + uint32_t *cigar = bam_get_cigar(b); + uint32_t end = bam_endpos(b); if (end >= tmp->start) { @@ -709,7 +757,6 @@ int readcore(struct depnode *header, bam1_t const *b, cntstat_t state) { if (isZero(tmp->len)) depnode_init(tmp); - tmp->cnts[pos - tmp->start]++; } for (i = 0; i < c->n_cigar; ++i) { @@ -724,8 +771,16 @@ int readcore(struct depnode *header, bam1_t const *b, cntstat_t state) continue; for (j = 0; j < l; ++j) { + /* once we have walked past the last region of this chromosome, + * the remaining bases of the CIGAR op (and later ops) cannot + * hit any region; just advance pos to stay consistent. */ + if (tmp == NULL) + { + pos += (l - j); + break; + } if (pos >= tmp->start) - match_pos(tmp, pos, tmp_state); + match_pos(&tmp, pos, tmp_state); pos++; } } @@ -806,12 +861,32 @@ int write_buffer_bgzf(kstring_t *str, BGZF *fp) return 0; } +/* Append one depth.tsv.gz row: "name\tpos\traw\trmdup\tcov\n". + * Uses kputw/kputuw instead of ksprintf to avoid the per-call vsnprintf + * overhead on the hot path (one call per target position). Output is + * byte-identical to ksprintf("%s\t%d\t%u\t%u\t%u\n", name,pos,raw,rmdup,cov). */ +static inline void kput_depth_line(kstring_t *s, const char *name, int name_len, + int pos, unsigned raw, unsigned rmdup, unsigned cov) +{ + kputsn(name, name_len, s); + kputc('\t', s); + kputw(pos, s); + kputc('\t', s); + kputuw(raw, s); + kputc('\t', s); + kputuw(rmdup, s); + kputc('\t', s); + kputuw(cov, s); + kputc('\n', s); +} + int stat_each_region(loopbams_parameters_t *para, aux_t *a) { struct depnode *node = para->tgt_node; if (isNull(node)) return 0; int j; + int name_len = para->name ? (int)strlen(para->name) : 0; float avg, med, cov1, cov2; uint32_t lst_start = 0; // uncover region start uint32_t lst_stop = 0; // uncover region stop @@ -824,8 +899,8 @@ int stat_each_region(loopbams_parameters_t *para, aux_t *a) cov2 = coverage_cal(node->covdep, node->len); for (j = 0; j < node->len; ++j) { - ksprintf(para->pdepths, "%s\t%d\t%u\t%u\t%u\n", para->name, node->start + j, node->vals[j], - node->rmdupdep[j], node->covdep[j]); + kput_depth_line(para->pdepths, para->name, name_len, (int)(node->start + j), + node->vals[j], node->rmdupdep[j], node->covdep[j]); // count_increase will alloc memory space automatically // use covdep to calculate coverage and averge depth count_increase(para->depvals_of_chr, node->covdep[j], uint32_t); @@ -855,27 +930,44 @@ int stat_each_region(loopbams_parameters_t *para, aux_t *a) { push_bedreg(para->ucreg, lst_start, lst_stop); } - write_buffer_bgzf(para->pdepths, para->fdep); + // 优化：只有当缓冲区超过阈值时才写入 + if (para->pdepths->l > WRITE_BUFFER_SIZE) + { + write_buffer_bgzf(para->pdepths, para->fdep); + } } else { + /* node never allocated (no read touched it): every position is zero + * depth. node->len == 0 here, so use the real span start..stop. This + * keeps depth.tsv.gz and the depth histograms consistent with regions + * that did receive reads. */ + unsigned length = node->stop - node->start + 1; avg = med = cov1 = cov2 = 0.0; - for (j = 0; j < node->len; ++j) + for (j = 0; j < (int)length; ++j) + { + kput_depth_line(para->pdepths, para->name, name_len, (int)(node->start + j), 0, 0, 0); + } + // 优化：批量写入而不是每行都写 + if (para->pdepths->l > WRITE_BUFFER_SIZE) { - ksprintf(para->pdepths, "%s\t%d\t0\t0\t0\n", para->name, node->start + j); write_buffer_bgzf(para->pdepths, para->fdep); } push_bedreg(para->ucreg, node->start, node->stop); // store uncover region - count_increaseN(para->depvals_of_chr, 0, node->len, uint32_t); - count_increaseN(a->c_dep, 0, node->len, uint32_t); + count_increaseN(para->depvals_of_chr, 0, length, uint32_t); + count_increaseN(a->c_dep, 0, length, uint32_t); + count_increaseN(a->c_rmdupdep, 0, length, uint32_t); } // ksprintf(para->pdepths,"\n"); count_increase(a->c_reg, (int)avg, uint32_t); ksprintf(para->rcov, "%s\t%u\t%u\t%.2f\t%.1f\t%.2f\t%.2f\n", para->name, node->start - 1, node->stop, avg, med, cov1, cov2); - // if (para->rcov->l > WINDOW_SIZE) - write_buffer_bgzf(para->rcov, para->freg); + // 优化：只有当缓冲区超过阈值时才写入 + if (para->rcov->l > WRITE_BUFFER_SIZE) + { + write_buffer_bgzf(para->rcov, para->freg); + } return 0; } @@ -907,11 +999,14 @@ int check_reachable_regions(loopbams_parameters_t *para, aux_t *a) ucreg_tmp = (bedreglist_t *)needmem(sizeof(bedreglist_t)); kh_val(h_uncov, l) = *ucreg_tmp; para->ucreg = &kh_val(h_uncov, l); + mustfree(ucreg_tmp); // data already copied into the hash table while (para->tgt_node) { - int length = para->tgt_node->stop - para->tgt_node->start + 1; - count_increaseN(a->c_dep, 0, length, uint32_t); + // stat_each_region now emits depth lines and feeds c_dep / + // c_rmdupdep / depvals_of_chr even for never-allocated nodes, + // so no manual count_increaseN is needed here (it would + // double-count the first, already-allocated node). stat_each_region(para, a); del_node(para->tgt_node); // no need allocate memory for these nodes } @@ -952,7 +1047,7 @@ void write_unover_file() int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) { // get the chromosome name from header - bam_header_t *h = a->h; + sam_hdr_t *h = a->h; loopbams_parameters_t *para = init_loopbams_parameters(); ksprintf(para->pdepths, "#Chr\tPos\tRaw Depth\tRmdup depth\tCover depth\n"); ksprintf(para->rcov, "#Chr\tStart\tStop\tAvg depth\tMedian\tCoverage\tCoverage(FIX)\n"); @@ -961,25 +1056,25 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) int i; for (i = 0; i < a->ndata; ++i) { - bamFile dat = a->data[i]; + htsFile *dat = a->data[i]; bool goto_next_chromosome = FALSE; int ret; cntstat_t state; // main loop bam1_t *b; - b = (bam1_t *)needmem(sizeof(bam1_t)); + b = bam_init1(); while (1) { state = CMATCH; - ret = bam_read1(dat, b); + ret = sam_read1(dat, h, b); if (ret == -1) { break; // normal end } - if (ret == -2) + if (ret < -1) { - errabort("%d bam file is truncated!\n", i + 1); + errabort("%d bam/cram file is truncated or corrupted!\n", i + 1); } bam1_core_t *c = &b->core; // skip secondary and supplement alignment first @@ -1031,7 +1126,17 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) // if tid != c->tid, completed the stat of the last chromosome // init the new chromosome node and clean the memory - // FIXME: need multi thread to improve it or NOT? + // + // 多线程优化说明： + // 当前实现是单线程顺序处理，可以考虑以下多线程方案： + // 1. 按染色体并行：每个线程处理一个染色体的深度统计 + // - 优点：实现简单，线程间无数据竞争 + // - 缺点：需要预先读取整个 BAM 文件或使用索引 + // 2. 生产者-消费者模式：一个线程读取 BAM，多个线程统计 + // - 优点：可以流式处理 + // - 缺点：需要线程安全的队列和同步机制 + // 目前保持单线程实现，因为 I/O 通常是瓶颈 + // if (para->tid != c->tid) { goto_next_chromosome = FALSE; // clean the flag of skiping @@ -1052,13 +1157,14 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) stat_flk_depcnt(para, a); while (para->flk_node) { - count_increaseN(a->c_flkdep, 0, para->flk_node->len, uint32_t); + int flen = para->flk_node->stop - para->flk_node->start + 1; + count_increaseN(a->c_flkdep, 0, flen, uint32_t); del_node(para->flk_node); } } para->tid = c->tid; - para->name = h->target_name[c->tid]; + para->name = (char *)sam_hdr_tid2name(h, c->tid); // impossible in normal pratices, only happens in the different bam // headers @@ -1093,15 +1199,20 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) para->flk_node = bed_depnode_list(para->flk); para->tar->flag = para->flk->flag = 1; - /* the next part is init uncover region hash*/ + /* 初始化 uncover region hash + * 注意：这里不能直接使用 ucreg_tmp 指针，因为 kh_val 返回的是值的拷贝 + * 我们需要获取 hash 表中实际存储位置的地址 + * 这样后续对 para->ucreg 的修改才能反映到 hash 表中 + */ k = kh_put(reg, h_uncov, strdup(para->name), &ret); bedreglist_t *ucreg_tmp; ucreg_tmp = (bedreglist_t *)needmem(sizeof(bedreglist_t)); + memset(ucreg_tmp, 0, sizeof(bedreglist_t)); // 确保初始化为零 kh_val(h_uncov, k) = *ucreg_tmp; - para->ucreg = &kh_val(h_uncov, k); // FIXME: I don't understand why - // couldn't use ucreg_tmp directly - /* finish init */ + mustfree(ucreg_tmp); // 释放临时变量，数据已拷贝到 hash 表 + para->ucreg = &kh_val(h_uncov, k); // 获取 hash 表中的实际地址 + /* finish init */ } while (para->flk_node && para->flk_node->stop < para->lstpos + 1) { @@ -1121,14 +1232,14 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) if (para->tgt_node && readcore(para->tgt_node, b, state)) { if (export_target_bam) - bam_write1(bamoutfp, b); + sam_write1(bamoutfp, a->h, b); fs->n_tgt++; } // endcore: } bam_destroy1(b); - bgzf_close(a->data[i]); + hts_close(a->data[i]); } while (para->tgt_node) { @@ -1137,7 +1248,8 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) } while (para->flk_node) { - count_increaseN(a->c_flkdep, 0, para->flk_node->len, uint32_t); + int flen = para->flk_node->stop - para->flk_node->start + 1; + count_increaseN(a->c_flkdep, 0, flen, uint32_t); del_node(para->flk_node); // no need allocate memory for these nodes } check_reachable_regions(para, a); @@ -1155,6 +1267,16 @@ struct regcov // MAX_CUTOFFS uint64_t cnt_array[10]; float cov_array[10]; + + // 新增：rmdup depth 覆盖度统计 + uint64_t cnt_rmdup, cnt4_rmdup, cnt10_rmdup, cnt30_rmdup, cnt100_rmdup; + float cov_rmdup, cov4_rmdup, cov10_rmdup, cov30_rmdup, cov100_rmdup; + uint64_t cnt_array_rmdup[10]; + float cov_array_rmdup[10]; + + // 新增：基于比例的覆盖度统计 + uint64_t cnt_ratio[10]; + float cov_ratio[10]; }; struct regcov *regcov_init() @@ -1225,7 +1347,9 @@ uint64_t cntcov_cal2(struct opt_aux *f, struct regcov *cov, count32_t *cnt, uint { (*data) += cnt->a[i] * i; } - uint64_t avg = (*data) / tgt_len; + // average depth (guard tgt_len==0 for zero-length BEDs; integer div-by-zero + // would crash the process) + uint64_t avg = tgt_len ? (*data) / tgt_len : 0; uint64_t avg_02 = 0.2 * avg; uint64_t avg_05 = 0.5 * avg; for (i = 0; i < cnt->m; ++i) @@ -1280,22 +1404,152 @@ uint64_t cntcov_cal2(struct opt_aux *f, struct regcov *cov, count32_t *cnt, uint return rawcnt; } -// FIXME: need improve soon!!! +// 计算 rmdup depth 的覆盖度统计 +uint64_t cntcov_cal_rmdup(struct opt_aux *f, struct regcov *cov, count32_t *cnt, uint64_t *data, uint64_t tgt_len) +{ + uint64_t rawcnt = 0; + int i; + *data = 0; + + // 计算总数据量和平均深度 + for (i = 0; i < cnt->m; ++i) + { + (*data) += cnt->a[i] * i; + rawcnt += cnt->a[i]; + } + + if (rawcnt == 0) + return 0; + + // 计算各阈值下的覆盖度 + uint64_t below_100 = 0, below_30 = 0, below_10 = 0, below_4 = 0; + + for (i = 0; i < cnt->m; ++i) + { + if (i < 100) + below_100 += cnt->a[i]; + if (i < 30) + below_30 += cnt->a[i]; + if (i < 10) + below_10 += cnt->a[i]; + if (i < 4) + below_4 += cnt->a[i]; + if (f->cutoff) + { + for (int x = 0; x < f->num_cutoffs; x++) + { + if (i < f->cutoffs[x]) + cov->cnt_array_rmdup[x] += cnt->a[i]; + } + } + } + + // 计算 rmdup 覆盖度 + cov->cnt_rmdup = rawcnt - (uint64_t)cnt->a[0]; + cov->cnt4_rmdup = rawcnt - below_4; + cov->cnt10_rmdup = rawcnt - below_10; + cov->cnt30_rmdup = rawcnt - below_30; + cov->cnt100_rmdup = rawcnt - below_100; + + cov->cov_rmdup = (float)cov->cnt_rmdup / rawcnt * 100; + cov->cov4_rmdup = (float)cov->cnt4_rmdup / rawcnt * 100; + cov->cov10_rmdup = (float)cov->cnt10_rmdup / rawcnt * 100; + cov->cov30_rmdup = (float)cov->cnt30_rmdup / rawcnt * 100; + cov->cov100_rmdup = (float)cov->cnt100_rmdup / rawcnt * 100; + + if (f->cutoff) + { + for (int x = 0; x < f->num_cutoffs; x++) + { + cov->cnt_array_rmdup[x] = rawcnt - cov->cnt_array_rmdup[x]; + cov->cov_array_rmdup[x] = (float)cov->cnt_array_rmdup[x] / rawcnt * 100; + } + } + + return rawcnt; +} + +// 计算基于比例的覆盖度统计 +void cntcov_cal_ratio(struct opt_aux *f, struct regcov *cov, count32_t *cnt, uint64_t tgt_len) +{ + if (!f->depth_ratio || f->num_ratios == 0) + return; + + // 计算总数据量 + uint64_t total_data = 0; + uint64_t rawcnt = 0; + int i; + + for (i = 0; i < cnt->m; ++i) + { + total_data += cnt->a[i] * i; + rawcnt += cnt->a[i]; + } + + if (rawcnt == 0 || tgt_len == 0) + return; + + // 计算平均深度 + float avg_depth = (float)total_data / tgt_len; + + // 对每个比例计算覆盖度 + for (int r = 0; r < f->num_ratios; r++) + { + uint64_t threshold = (uint64_t)(f->ratios[r] * avg_depth); + uint64_t below_threshold = 0; + + for (i = 0; i < cnt->m && i <= threshold; ++i) + { + below_threshold += cnt->a[i]; + } + + cov->cnt_ratio[r] = rawcnt - below_threshold; + cov->cov_ratio[r] = (float)cov->cnt_ratio[r] / rawcnt * 100; + } +} + +// 计算深度分布的中位数 +// 使用累积分布查找，时间复杂度 O(n) float median_cnt(count32_t *cnt) { + if (cnt->m == 0) + return 0; + int i; uint64_t sum = 0; + + // 计算总数 for (i = 0; i < cnt->m; ++i) sum += (uint64_t)cnt->a[i]; - uint64_t med = sum / 2; + + if (sum == 0) + return 0; + + // 1-indexed middle positions: for odd sum, (sum+1)/2 is the single middle; + // for even sum, sum/2 and sum/2+1 are the two middles. The depth at + // position p is the smallest i with cumulative count >= p. + uint64_t lower = (sum + 1) / 2; // lower-middle (== upper when sum is odd) + uint64_t upper = sum / 2 + 1; // upper-middle uint64_t num = 0; + int v1 = -1, v2 = -1; for (i = 0; i < cnt->m; ++i) { num += (uint64_t)cnt->a[i]; - if (num >= med) - return (float)i; + if (v1 < 0 && num >= lower) + v1 = i; + if (v2 < 0 && num >= upper) + { + v2 = i; + break; + } } - return 0; + if (v1 < 0) + v1 = cnt->m - 1; // fallback, should not happen + if (v2 < 0) + v2 = v1; + if (sum & 1) + return (float)v1; + return (float)(v1 + v2) / 2.0f; } float average_cnt(count32_t *cnt) @@ -1308,7 +1562,273 @@ float average_cnt(count32_t *cnt) sum += (uint64_t)cnt->a[i] * i; num += (uint64_t)cnt->a[i]; } - return (float)sum / num; + return num ? (float)sum / num : 0; +} + +// JSON 格式化输出 - 缩进级别跟踪 +static int json_indent_level = 0; +static const char *JSON_INDENT = " "; // 2空格缩进 + +static void json_newline_indent(kstring_t *s) +{ + kputc('\n', s); + for (int i = 0; i < json_indent_level; i++) + kputs(JSON_INDENT, s); +} + +static void json_start_object(kstring_t *s) +{ + kputc('{', s); + json_indent_level++; +} + +static void json_end_object(kstring_t *s) +{ + // Remove trailing comma if present + if (s->l > 0 && s->s[s->l - 1] == ',') + s->l--; + json_indent_level--; + json_newline_indent(s); + kputc('}', s); +} + +static void json_start_array(kstring_t *s) +{ + kputc('[', s); +} + +static void json_end_array(kstring_t *s) +{ + if (s->l > 0 && s->s[s->l - 1] == ',') + s->l--; + kputc(']', s); +} + +static void json_key(kstring_t *s, const char *key) +{ + json_newline_indent(s); + ksprintf(s, "\"%s\": ", key); +} + +static void json_string(kstring_t *s, const char *key, const char *value) +{ + json_newline_indent(s); + ksprintf(s, "\"%s\": \"%s\",", key, value); +} + +static void json_int(kstring_t *s, const char *key, int64_t value) +{ + json_newline_indent(s); + ksprintf(s, "\"%s\": %" PRId64 ",", key, value); +} + +static void json_uint(kstring_t *s, const char *key, uint64_t value) +{ + json_newline_indent(s); + ksprintf(s, "\"%s\": %" PRIu64 ",", key, value); +} + +static void json_float(kstring_t *s, const char *key, float value) +{ + json_newline_indent(s); + ksprintf(s, "\"%s\": %.2f,", key, value); +} + +/* percentage guarded against a zero denominator (empty BAM / all-unmapped / + * zero-length target). Returns 0.0 instead of NaN/Inf. */ +static inline float safe_pct(uint64_t num, uint64_t den) +{ + return den ? (float)num / (float)den * 100 : 0.0f; +} + +/* Generate JSON format coverage report */ +int print_report_json(struct opt_aux *f, aux_t *a, bamflag_t *fs, + struct regcov *tarcov, struct regcov *flkcov, + struct regcov *regcov, float iavg, uint64_t imed) +{ + if (outdir) chdir(outdir); + FILE *fjson = open_wfile("coverage.report.json"); + if (!fjson) { + warnings("Failed to create coverage.report.json"); + return -1; + } + // 设置更大的文件缓冲区 + setvbuf(fjson, NULL, _IOFBF, WRITE_BUFFER_SIZE); + + // 重置缩进级别 + json_indent_level = 0; + + kstring_t json = {0, 0, NULL}; + + json_start_object(&json); + + // Version and files + json_string(&json, "version", Version); + json_key(&json, "files"); + json_start_array(&json); + for (int i = 0; i < f->nfiles; ++i) { + ksprintf(&json, "\"%s\",", f->inputs[i]); + } + json_end_array(&json); + kputc(',', &json); + + // Total section + json_key(&json, "total"); + json_start_object(&json); + json_uint(&json, "raw_reads", fs->n_reads); + json_uint(&json, "qc_fail_reads", fs->n_qcfail); + json_float(&json, "raw_data_mb", (float)fs->n_data / 1e6); + json_uint(&json, "paired_reads", fs->n_pair_all); + json_uint(&json, "mapped_reads", fs->n_mapped); + json_float(&json, "mapped_reads_fraction", safe_pct(fs->n_mapped, fs->n_reads)); + json_float(&json, "mapped_data_mb", (float)fs->n_mdata / 1e6); + json_float(&json, "mapped_data_fraction", safe_pct(fs->n_mdata, fs->n_data)); + json_uint(&json, "properly_paired", fs->n_pair_good); + json_float(&json, "properly_paired_fraction", safe_pct(fs->n_pair_good, fs->n_reads)); + json_uint(&json, "read_mate_paired", fs->n_pair_map); + json_float(&json, "read_mate_paired_fraction", safe_pct(fs->n_pair_map, fs->n_reads)); + json_uint(&json, "singletons", fs->n_sgltn); + json_uint(&json, "diff_chr", fs->n_diffchr); + json_uint(&json, "read1", fs->n_read1); + json_uint(&json, "read2", fs->n_read2); + json_uint(&json, "read1_rmdup", fs->n_rmdup1); + json_uint(&json, "read2_rmdup", fs->n_rmdup2); + json_uint(&json, "forward_strand", fs->n_pstrand); + json_uint(&json, "backward_strand", fs->n_mstrand); + json_uint(&json, "pcr_duplicates", fs->n_dup); + json_float(&json, "pcr_duplicates_fraction", safe_pct(fs->n_dup, fs->n_mapped)); + json_int(&json, "mapq_cutoff", f->mapQ_lim); + json_uint(&json, "mapq_reads", fs->n_qual); + json_float(&json, "mapq_reads_fraction_all", safe_pct(fs->n_qual, fs->n_reads)); + json_float(&json, "mapq_reads_fraction_mapped", safe_pct(fs->n_qual, fs->n_mapped)); + json_end_object(&json); + kputc(',', &json); + + // Insert size section + json_key(&json, "insert_size"); + json_start_object(&json); + json_float(&json, "average", iavg); + json_uint(&json, "median", imed); + json_end_object(&json); + kputc(',', &json); + + // Target section + json_key(&json, "target"); + json_start_object(&json); + json_uint(&json, "target_reads", fs->n_tgt); + json_float(&json, "target_reads_fraction_all", safe_pct(fs->n_tgt, fs->n_reads)); + json_float(&json, "target_reads_fraction_mapped", safe_pct(fs->n_tgt, fs->n_mapped)); + json_float(&json, "target_data_mb", (float)fs->n_tdata / 1e6); + json_float(&json, "target_data_rmdup_mb", (float)fs->n_trmdat / 1e6); + json_float(&json, "target_data_fraction_all", safe_pct(fs->n_tdata, fs->n_data)); + json_float(&json, "target_data_fraction_mapped", safe_pct(fs->n_tdata, fs->n_mdata)); + json_uint(&json, "region_length", a->tgt_len); + json_float(&json, "average_depth", (float)fs->n_tdata / (a->tgt_len ? a->tgt_len : 1)); + json_float(&json, "average_depth_rmdup", (float)fs->n_trmdat / (a->tgt_len ? a->tgt_len : 1)); + + // Coverage sub-object + json_key(&json, "coverage"); + json_start_object(&json); + json_float(&json, "gt_0x", tarcov->cov); + json_float(&json, "gte_4x", tarcov->cov4); + json_float(&json, "gte_10x", tarcov->cov10); + json_float(&json, "gte_30x", tarcov->cov30); + json_float(&json, "gte_100x", tarcov->cov100); + json_float(&json, "gt_0_2_avg", tarcov->cov02x); + json_float(&json, "gt_0_5_avg", tarcov->cov05x); + if (f->cutoff) { + json_key(&json, "custom"); + json_start_object(&json); + for (int x = 0; x < f->num_cutoffs; x++) { + char key[32]; + sprintf(key, "gte_%dx", f->cutoffs[x]); + json_float(&json, key, tarcov->cov_array[x]); + } + json_end_object(&json); + kputc(',', &json); + } + if (f->depth_ratio && f->num_ratios > 0) { + json_key(&json, "ratio_based"); + json_start_object(&json); + for (int r = 0; r < f->num_ratios; r++) { + char key[32]; + sprintf(key, "gt_%.1f_avg", f->ratios[r]); + json_float(&json, key, tarcov->cov_ratio[r]); + } + json_end_object(&json); + kputc(',', &json); + } + json_end_object(&json); + kputc(',', &json); + + // Coverage rmdup sub-object + json_key(&json, "coverage_rmdup"); + json_start_object(&json); + json_float(&json, "gt_0x", tarcov->cov_rmdup); + json_float(&json, "gte_4x", tarcov->cov4_rmdup); + json_float(&json, "gte_10x", tarcov->cov10_rmdup); + json_float(&json, "gte_30x", tarcov->cov30_rmdup); + json_float(&json, "gte_100x", tarcov->cov100_rmdup); + if (f->cutoff) { + json_key(&json, "custom"); + json_start_object(&json); + for (int x = 0; x < f->num_cutoffs; x++) { + char key[32]; + sprintf(key, "gte_%dx", f->cutoffs[x]); + json_float(&json, key, tarcov->cov_array_rmdup[x]); + } + json_end_object(&json); + kputc(',', &json); + } + json_end_object(&json); + kputc(',', &json); + + // Region coverage + json_uint(&json, "region_count", a->tgt_nreg); + json_key(&json, "region_coverage"); + json_start_object(&json); + json_uint(&json, "gt_0x_count", regcov->cnt); + json_float(&json, "gt_0x_fraction", regcov->cov); + json_float(&json, "gte_4x_fraction", regcov->cov4); + json_float(&json, "gte_10x_fraction", regcov->cov10); + json_float(&json, "gte_30x_fraction", regcov->cov30); + json_float(&json, "gte_100x_fraction", regcov->cov100); + json_end_object(&json); + + json_end_object(&json); // end target + kputc(',', &json); + + // Flank section + json_key(&json, "flank"); + json_start_object(&json); + json_int(&json, "flank_size", flank_reg); + json_uint(&json, "region_length", a->flk_len); + json_float(&json, "average_depth", (float)fs->n_fdata / (a->flk_len ? a->flk_len : 1)); + json_uint(&json, "flank_reads", fs->n_flk); + json_float(&json, "flank_reads_fraction_all", safe_pct(fs->n_flk, fs->n_reads)); + json_float(&json, "flank_reads_fraction_mapped", safe_pct(fs->n_flk, fs->n_mapped)); + json_float(&json, "flank_data_mb", (float)fs->n_fdata / 1e6); + json_float(&json, "flank_data_fraction_all", safe_pct(fs->n_fdata, fs->n_data)); + json_float(&json, "flank_data_fraction_mapped", safe_pct(fs->n_fdata, fs->n_mdata)); + json_key(&json, "coverage"); + json_start_object(&json); + json_float(&json, "gt_0x", flkcov->cov); + json_float(&json, "gte_4x", flkcov->cov4); + json_float(&json, "gte_10x", flkcov->cov10); + json_float(&json, "gte_30x", flkcov->cov30); + json_float(&json, "gte_100x", flkcov->cov100); + json_end_object(&json); + + json_end_object(&json); // end flank + + json_end_object(&json); // end root + + // Write to file + fprintf(fjson, "%s\n", json.s); + fclose(fjson); + free(json.s); + + return 0; } int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) @@ -1319,6 +1839,9 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) FILE *fdep; finsert = open_wfile("insertsize.plot"); fdep = open_wfile("depth_distribution.plot"); + // 设置更大的文件缓冲区以减少系统调用 + setvbuf(finsert, NULL, _IOFBF, WRITE_BUFFER_SIZE); + setvbuf(fdep, NULL, _IOFBF, WRITE_BUFFER_SIZE); struct regcov *tarcov = regcov_init(); struct regcov *flkcov = regcov_init(); @@ -1357,8 +1880,17 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) cntcov_cal(f, flkcov, a->c_flkdep, &fs->n_fdata); for (i = 0; i < a->c_rmdupdep->m; ++i) fs->n_trmdat += a->c_rmdupdep->a[i] * i; + + // 计算 rmdup depth 覆盖度统计 + uint64_t rmdup_data; + cntcov_cal_rmdup(f, tarcov, a->c_rmdupdep, &rmdup_data, a->tgt_len); + + // 计算基于比例的覆盖度统计 + cntcov_cal_ratio(f, tarcov, a->c_dep, a->tgt_len); FILE *fchrcov = open_wfile("chromosomes.report"); + // 设置更大的文件缓冲区 + setvbuf(fchrcov, NULL, _IOFBF, WRITE_BUFFER_SIZE); { fprintf(fchrcov, "%11s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s", "#Chromosome", "DATA(%)", "Avg depth", "Median", "Coverage%", "Cov 4x %", "Cov 10x %", "Cov 30x %", "Cov 100x %"); @@ -1419,6 +1951,8 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) } fclose(fchrcov); FILE *fc = open_wfile("coverage.report"); + // 设置更大的文件缓冲区以减少系统调用 + setvbuf(fc, NULL, _IOFBF, WRITE_BUFFER_SIZE); do { fprintf(fc, "## The file was created by %s\n", program_name); @@ -1432,15 +1966,15 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) fprintf(fc, "%60s\t%.2f\n", "[Total] Raw Data(Mb)", (float)fs->n_data / 1e6); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] Paired Reads", fs->n_pair_all); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] Mapped Reads", fs->n_mapped); - fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of Mapped Reads", (float)fs->n_mapped / fs->n_reads * 100); + fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of Mapped Reads", safe_pct(fs->n_mapped, fs->n_reads)); fprintf(fc, "%60s\t%.2f\n", "[Total] Mapped Data(Mb)", fs->n_mdata / 1e6); - fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of Mapped Data(Mb)", (float)fs->n_mdata / fs->n_data * 100); + fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of Mapped Data(Mb)", safe_pct(fs->n_mdata, fs->n_data)); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] Properly paired", fs->n_pair_good); fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of Properly paired", - (float)fs->n_pair_good / fs->n_reads * 100); + safe_pct(fs->n_pair_good, fs->n_reads)); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] Read and mate paired", fs->n_pair_map); fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of Read and mate paired", - (float)fs->n_pair_map / fs->n_reads * 100); + safe_pct(fs->n_pair_map, fs->n_reads)); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] Singletons", fs->n_sgltn); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] Read and mate map to diff chr", fs->n_diffchr); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] Read1", fs->n_read1); @@ -1451,31 +1985,31 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] backward strand reads", fs->n_mstrand); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] PCR duplicate reads", fs->n_dup); fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of PCR duplicate reads", - (float)fs->n_dup / fs->n_mapped * 100); // change n_reads to n_mapped, 2015/05/25 + safe_pct(fs->n_dup, fs->n_mapped)); // change n_reads to n_mapped, 2015/05/25 fprintf(fc, "%60s\t%d\n", "[Total] Map quality cutoff value", f->mapQ_lim); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Total] MapQuality above cutoff reads", fs->n_qual); fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of MapQ reads in all reads", - (float)fs->n_qual / fs->n_reads * 100); + safe_pct(fs->n_qual, fs->n_reads)); fprintf(fc, "%60s\t%.2f%%\n", "[Total] Fraction of MapQ reads in mapped reads", - (float)fs->n_qual / fs->n_mapped * 100); + safe_pct(fs->n_qual, fs->n_mapped)); // insert fprintf(fc, "%60s\t%.2f\n", "[Insert size] Average", iavg); fprintf(fc, "%60s\t%.ld\n", "[Insert size] Median", imed); // tgt fprintf(fc, "%60s\t%" PRIu64 "\n", "[Target] Target Reads", fs->n_tgt); fprintf(fc, "%60s\t%.2f%%\n", "[Target] Fraction of Target Reads in all reads", - (float)fs->n_tgt / fs->n_reads * 100); + safe_pct(fs->n_tgt, fs->n_reads)); fprintf(fc, "%60s\t%.2f%%\n", "[Target] Fraction of Target Reads in mapped reads", - (float)fs->n_tgt / fs->n_mapped * 100); + safe_pct(fs->n_tgt, fs->n_mapped)); fprintf(fc, "%60s\t%.2f\n", "[Target] Target Data(Mb)", (float)fs->n_tdata / 1e6); fprintf(fc, "%60s\t%.2f\n", "[Target] Target Data Rmdup(Mb)", (float)fs->n_trmdat / 1e6); fprintf(fc, "%60s\t%.2f%%\n", "[Target] Fraction of Target Data in all data", - (float)fs->n_tdata / fs->n_data * 100); + safe_pct(fs->n_tdata, fs->n_data)); fprintf(fc, "%60s\t%.2f%%\n", "[Target] Fraction of Target Data in mapped data", - (float)fs->n_tdata / fs->n_mdata * 100); + safe_pct(fs->n_tdata, fs->n_mdata)); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Target] Len of region", a->tgt_len); - fprintf(fc, "%60s\t%.2f\n", "[Target] Average depth", (float)fs->n_tdata / a->tgt_len); - fprintf(fc, "%60s\t%.2f\n", "[Target] Average depth(rmdup)", (float)fs->n_trmdat / a->tgt_len); + fprintf(fc, "%60s\t%.2f\n", "[Target] Average depth", (float)fs->n_tdata / (a->tgt_len ? a->tgt_len : 1)); + fprintf(fc, "%60s\t%.2f\n", "[Target] Average depth(rmdup)", (float)fs->n_trmdat / (a->tgt_len ? a->tgt_len : 1)); fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage (>0.2*(Average depth)x)", tarcov->cov02x); fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage (>0.5*(Average depth)x)", tarcov->cov05x); fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage (>0x)", tarcov->cov); @@ -1492,6 +2026,31 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) fprintf(fc, "%60s\t%.2f%%\n", titles, tarcov->cov_array[x]); } } + // rmdup coverage statistics + fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage(rmdup) (>0x)", tarcov->cov_rmdup); + fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage(rmdup) (>=4x)", tarcov->cov4_rmdup); + fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage(rmdup) (>=10x)", tarcov->cov10_rmdup); + fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage(rmdup) (>=30x)", tarcov->cov30_rmdup); + fprintf(fc, "%60s\t%.2f%%\n", "[Target] Coverage(rmdup) (>=100x)", tarcov->cov100_rmdup); + if (f->cutoff) + { + char titles[60]; + for (int x = 0; x < f->num_cutoffs; x++) + { + sprintf(titles, "[Target] Coverage(rmdup) (>=%ux)", f->cutoffs[x]); + fprintf(fc, "%60s\t%.2f%%\n", titles, tarcov->cov_array_rmdup[x]); + } + } + // ratio-based coverage statistics + if (f->depth_ratio && f->num_ratios > 0) + { + char titles[60]; + for (int r = 0; r < f->num_ratios; r++) + { + sprintf(titles, "[Target] Coverage (>%.1f*Avg)", f->ratios[r]); + fprintf(fc, "%60s\t%.2f%%\n", titles, tarcov->cov_ratio[r]); + } + } // tgt regions fprintf(fc, "%60s\t%u\n", "[Target] Target Region Count", a->tgt_nreg); fprintf(fc, "%60s\t%" PRIu64 "\n", "[Target] Region covered > 0x", regcov->cnt); @@ -1514,17 +2073,17 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) // flk fprintf(fc, "%60s\t%u\n", "[flank] flank size", flank_reg); fprintf(fc, "%60s\t%" PRIu64 "\n", "[flank] Len of region (not include target region)", a->flk_len); - fprintf(fc, "%60s\t%.2f\n", "[flank] Average depth", (float)fs->n_fdata / a->flk_len); + fprintf(fc, "%60s\t%.2f\n", "[flank] Average depth", (float)fs->n_fdata / (a->flk_len ? a->flk_len : 1)); fprintf(fc, "%60s\t%" PRIu64 "\n", "[flank] flank Reads", fs->n_flk); fprintf(fc, "%60s\t%.2f%%\n", "[flank] Fraction of flank Reads in all reads", - (float)fs->n_flk / fs->n_reads * 100); + safe_pct(fs->n_flk, fs->n_reads)); fprintf(fc, "%60s\t%.2f%%\n", "[flank] Fraction of flank Reads in mapped reads", - (float)fs->n_flk / fs->n_mapped * 100); + safe_pct(fs->n_flk, fs->n_mapped)); fprintf(fc, "%60s\t%.2f\n", "[flank] flank Data(Mb)", (float)fs->n_fdata / 1e6); fprintf(fc, "%60s\t%.2f%%\n", "[flank] Fraction of flank Data in all data", - (float)fs->n_fdata / fs->n_data * 100); + safe_pct(fs->n_fdata, fs->n_data)); fprintf(fc, "%60s\t%.2f%%\n", "[flank] Fraction of flank Data in mapped data", - (float)fs->n_fdata / fs->n_mdata * 100); + safe_pct(fs->n_fdata, fs->n_mdata)); fprintf(fc, "%60s\t%.2f%%\n", "[flank] Coverage (>0x)", flkcov->cov); fprintf(fc, "%60s\t%.2f%%\n", "[flank] Coverage (>=4x)", flkcov->cov4); fprintf(fc, "%60s\t%.2f%%\n", "[flank] Coverage (>=10x)", flkcov->cov10); @@ -1542,6 +2101,9 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) } while (0); + // 生成 JSON 格式报告 + print_report_json(f, a, fs, tarcov, flkcov, regcov, iavg, imed); + mustfree(tarcov); mustfree(regcov); mustfree(flkcov); @@ -1556,6 +2118,9 @@ enum INSERTSIZE, UNCOVER, BAMOUT, + DEPTHRATIO, + REFERENCE, + THREADS, HELP }; @@ -1568,6 +2133,9 @@ static struct option const long_opts[] = {{"outdir", required_argument, NULL, 'o {"mapthres", required_argument, NULL, 'q'}, {"uncover", required_argument, NULL, UNCOVER}, {"bamout", required_argument, NULL, BAMOUT}, + {"depthratio", required_argument, NULL, DEPTHRATIO}, + {"reference", required_argument, NULL, 'T'}, + {"threads", required_argument, NULL, THREADS}, //{"rmdup", no_argument, NULL, 'd'}, {"help", no_argument, NULL, 'h'}, {"version", no_argument, NULL, 'v'}}; @@ -1577,14 +2145,14 @@ int show_version() printf("%s\n", Version); return 1; } -int bamdst(int argc, char *argv[]) +int xamdst(int argc, char *argv[]) { int n, i; char *probe = 0; // struct opt_aux opt = {.inputs = NULL, .isize_lim = 2000, .mapQ_lim = 20}; struct opt_aux opt = init_opt_aux(); - while ((n = getopt_long(argc, argv, "o:p:f:q:l:h1v", long_opts, NULL)) >= 0) + while ((n = getopt_long(argc, argv, "o:p:f:q:l:T:h1v", long_opts, NULL)) >= 0) { switch (n) { @@ -1632,9 +2200,44 @@ int bamdst(int argc, char *argv[]) case BAMOUT: export_target_bam = strdup(optarg); break; + case DEPTHRATIO: + // 解析用户指定的深度比例 + { + char *ratio_token = strtok(optarg, ","); + opt.num_ratios = 0; // 重置为用户指定的值 + while (ratio_token != NULL) + { + if (opt.num_ratios >= opt.max_ratios) + { + fprintf(stderr, "Too many depth ratios specified (max %d)\n", opt.max_ratios); + exit(EXIT_FAILURE); + } + float ratio = atof(ratio_token); + if (ratio <= 0.0f || ratio > 1.0f) + { + fprintf(stderr, "Invalid depth ratio: %s (must be between 0 and 1)\n", ratio_token); + exit(EXIT_FAILURE); + } + opt.ratios[opt.num_ratios++] = ratio; + ratio_token = strtok(NULL, ","); + } + } + break; case 'q': opt.mapQ_lim = atoi(optarg); break; + case 'T': + // 参考基因组路径（用于 CRAM 文件） + opt.reference = strdup(optarg); + break; + case THREADS: + opt.nthreads = atoi(optarg); + if (opt.nthreads < 0) + { + fprintf(stderr, "Invalid thread count: %s\n", optarg); + exit(EXIT_FAILURE); + } + break; case 'h': usage(1); break; @@ -1650,7 +2253,7 @@ int bamdst(int argc, char *argv[]) } } if (isNull(outdir) || isNull(probe)) - usage(0); + usage(1); if (export_target_bam && check_filename_isbam(export_target_bam)) { fprintf(stderr, "--bamout must be a bam file: %s", export_target_bam); @@ -1659,67 +2262,136 @@ int bamdst(int argc, char *argv[]) n = argc - optind; mkdirp(outdir, 0755); - // capable of deals with severl bam files + // capable of deals with several bam/cram files aux_t *aux; aux = aux_init(); if (isZero(n)) { - aux->data = (bamFile *)needmem(sizeof(bamFile)); - aux->data[0] = bgzf_dopen(fileno(stdin), "r"); - aux->h = bam_header_read(aux->data[0]); + aux->data = (htsFile **)needmem(sizeof(htsFile *)); + aux->data[0] = hts_open("-", "r"); + if (aux->data[0] == NULL) + errabort("Failed to open stdin for reading"); + // Enable multi-threading for stdin reading + if (opt.nthreads > 0) + { + if (hts_set_threads(aux->data[0], opt.nthreads) != 0) + { + fprintf(stderr, "Warning: Failed to set %d threads for stdin\n", opt.nthreads); + } + } + // Set reference for CRAM if provided + if (opt.reference) + { + if (hts_set_fai_filename(aux->data[0], opt.reference) != 0) + { + errabort("Failed to set reference file: %s", opt.reference); + } + } + aux->h = sam_hdr_read(aux->data[0]); + if (aux->h == NULL) + errabort("Failed to read header from stdin"); aux->ndata = 1; opt.nfiles = 0; } else { - aux->data = (bamFile *)needmem(n * sizeof(bamFile)); + aux->data = (htsFile **)needmem(n * sizeof(htsFile *)); opt.nfiles = n; opt.inputs = (char **)needmem(n * sizeof(char *)); for (i = 0; i < n; ++i) { - bam_header_t *h_tmp; - // h_tmp = calloc(1, sizeof(bam_header_t)); + sam_hdr_t *h_tmp; if (STREQ(argv[optind + i], "-")) { - aux->data[i] = bgzf_dopen(fileno(stdin), "r"); + aux->data[i] = hts_open("-", "r"); stdin_lock = 1; } else { - aux->data[i] = bgzf_open(argv[optind + i], "r"); + aux->data[i] = hts_open(argv[optind + i], "r"); } if (aux->data[i] == NULL) errabort("%s: %s", argv[optind + i], strerror(errno)); - h_tmp = bam_header_read(aux->data[i]); + + // Enable multi-threading for BAM/CRAM reading + if (opt.nthreads > 0) + { + if (hts_set_threads(aux->data[i], opt.nthreads) != 0) + { + fprintf(stderr, "Warning: Failed to set %d threads for %s\n", opt.nthreads, argv[optind + i]); + } + } + + // Check if CRAM and reference is needed + const htsFormat *fmt = hts_get_format(aux->data[i]); + if (fmt->format == cram) + { + if (opt.reference == NULL) + { + errabort("CRAM file requires a reference genome. Use -T/--reference to specify."); + } + if (hts_set_fai_filename(aux->data[i], opt.reference) != 0) + { + errabort("Failed to set reference file: %s", opt.reference); + } + } + else if (opt.reference) + { + // Also set reference for BAM if provided (useful for some operations) + hts_set_fai_filename(aux->data[i], opt.reference); + } + + h_tmp = sam_hdr_read(aux->data[i]); + if (h_tmp == NULL) + errabort("Failed to read header from %s", argv[optind + i]); if (i == 0) aux->h = h_tmp; else - bam_header_destroy(h_tmp); + sam_hdr_destroy(h_tmp); opt.inputs[i] = strdup(argv[optind + i]); } aux->ndata = n; } - // FIXME: accpet more than one bam files! + // 多 BAM/CRAM 文件支持说明： + // 当前实现已支持多个 BAM/CRAM 文件输入，它们会被顺序处理并合并统计结果 + // 注意：所有文件必须使用相同的参考基因组和排序方式 + // 第一个文件的 header 会被用于输出 if (export_target_bam) { - bamoutfp = bam_open(export_target_bam, "w"); + bamoutfp = hts_open(export_target_bam, "wb"); if (bamoutfp == NULL) errabort("%s : %s", export_target_bam, strerror(errno)); - bam_header_write(bamoutfp, aux->h); + // Enable multi-threading for BAM writing + if (opt.nthreads > 0) + { + if (hts_set_threads(bamoutfp, opt.nthreads) != 0) + { + fprintf(stderr, "Warning: Failed to set %d threads for output BAM\n", opt.nthreads); + } + } + // Write BAM header using htslib compatible format + if (sam_hdr_write(bamoutfp, aux->h) < 0) + errabort("Failed to write header to %s", export_target_bam); } h_chrlength_init(); header2chrhash(aux->h); load_bed_init(probe, aux); chrhash_destroy(); freemem(probe); - if (aux->c_isize->n < opt.isize_lim) + // Pre-size the insert-size histogram so count_increase does not trigger a + // realloc (and a transient shrink) on the first large insert sizes. Must + // keep n in sync with the real allocation. + if (opt.isize_lim > 0 && aux->c_isize->n < (unsigned)opt.isize_lim) { - aux->c_isize->a = realloc(aux->c_isize->a, opt.isize_lim * sizeof(unsigned)); + unsigned *isize_tmp = realloc(aux->c_isize->a, (size_t)opt.isize_lim * sizeof(unsigned)); + if (isize_tmp == NULL) + errabort("out of memory while resizing insert-size histogram"); + aux->c_isize->a = isize_tmp; + aux->c_isize->n = opt.isize_lim; } for (i = 0; i < opt.isize_lim; ++i) aux->c_isize->a[i] = 0; - // aux->c_isize->m = opt.isize_lim; - aux->nchr = aux->h->n_targets; + aux->nchr = sam_hdr_nref(aux->h); struct bamflag fs = {}; load_bamfiles(&opt, aux, &fs); print_report(&opt, aux, &fs); @@ -1728,8 +2400,10 @@ int bamdst(int argc, char *argv[]) freemem(opt.inputs[i]); freemem(opt.inputs); if (export_target_bam) - bam_close(bamoutfp); + hts_close(bamoutfp); freemem(opt.cutoffs); + freemem(opt.ratios); + freemem(opt.reference); freeall: freemem(export_target_bam); freemem(outdir); @@ -1739,5 +2413,5 @@ int bamdst(int argc, char *argv[]) /* main */ int main(int argc, char *argv[]) { - return bamdst(argc, argv); + return xamdst(argc, argv); }