From e40dc823c3f38c6407e271b80698fa2e9c369c91 Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 24 Dec 2025 13:58:07 +0800 Subject: [PATCH 01/18] =?UTF-8?q?[update]=20=E6=96=B0=E5=A2=9Ejson?= =?UTF-8?q?=E7=BB=93=E6=9E=9C=E5=92=8C=E6=9B=B4=E5=A4=9A=E7=9A=84=E8=A6=86?= =?UTF-8?q?=E7=9B=96=E5=BA=A6=E7=BB=9F=E8=AE=A1=E6=A8=A1=E5=BC=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- README.md | 12 +- bamdst.c | 443 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 451 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 99efe9c..0930894 100644 --- a/README.md +++ b/README.md @@ -71,18 +71,22 @@ value for reasonal visual purpose. Default is 2000. set this cutoff value for calculate the bad covered region. Default is <5. ---use_rmdup (an invalid parament since v1.0.0 ) +--depthratio [ratios] -Use rmdup depth instead of cover depth to calculate the coverage of target regions and +specify depth ratios for coverage calculation relative to average depth. -so on. +Use comma-separated values like "0.1,0.2,0.5". Default is "0.2,0.5". + +This calculates coverage at positions with depth > (ratio × average_depth). ## OUTPUT FILES -Seven files will be created in the output direction. There are: +Eight files will be created in the output direction. There are: -**coverage.report** +-**coverage.report.json** (NEW: JSON format for programmatic access) + -**cumu.plot** -**insert.plot** diff --git a/bamdst.c b/bamdst.c index 285b15e..9450ecc 100644 --- a/bamdst.c +++ b/bamdst.c @@ -154,6 +154,11 @@ 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) }; // 一个函数来初始化 opt_aux 结构体 @@ -175,6 +180,21 @@ 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; + return opt; } @@ -458,6 +478,7 @@ Optional parameters:\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\ @@ -471,6 +492,7 @@ Optional parameters:\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\ @@ -478,6 +500,11 @@ Optional parameters:\n\ - depth.tsv.gz raw depth, rmdup depth, coverage depth of each position\n\ - uncover.bed the bad covered or uncovered region in the probe file\n\ \n\ +* New features in this version:\n\ + - Rmdup coverage statistics: coverage based on deduplicated depth\n\ + - Ratio-based coverage: coverage at ratios of average depth (e.g., 0.2x, 0.5x)\n\ + - JSON output: coverage.report.json for programmatic access\n\ +\n\ * About depth.tsv.gz:\n\ * There are five columns in this file, including chromosome, position, raw\n\ * depth, rmdep depth, coverage depth\n\ @@ -1155,6 +1182,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() @@ -1280,6 +1317,110 @@ uint64_t cntcov_cal2(struct opt_aux *f, struct regcov *cov, count32_t *cnt, uint return rawcnt; } +// 计算 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, below_0 = 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; + } +} + // FIXME: need improve soon!!! float median_cnt(count32_t *cnt) { @@ -1311,6 +1452,247 @@ float average_cnt(count32_t *cnt) return (float)sum / num; } +/* JSON string builder utilities */ +static void json_start_object(kstring_t *s) +{ + kputc('{', s); +} + +static void json_end_object(kstring_t *s) +{ + // Remove trailing comma if present + if (s->l > 0 && s->s[s->l - 1] == ',') + s->l--; + 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) +{ + ksprintf(s, "\"%s\":", key); +} + +static void json_string(kstring_t *s, const char *key, const char *value) +{ + ksprintf(s, "\"%s\":\"%s\",", key, value); +} + +static void json_int(kstring_t *s, const char *key, int64_t value) +{ + ksprintf(s, "\"%s\":%" PRId64 ",", key, value); +} + +static void json_uint(kstring_t *s, const char *key, uint64_t value) +{ + ksprintf(s, "\"%s\":%" PRIu64 ",", key, value); +} + +static void json_float(kstring_t *s, const char *key, float value) +{ + ksprintf(s, "\"%s\":%.2f,", key, value); +} + +static void json_float_array_value(kstring_t *s, float value) +{ + ksprintf(s, "%.2f,", value); +} + +/* 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; + } + + 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", (float)fs->n_mapped / fs->n_reads * 100); + json_float(&json, "mapped_data_mb", (float)fs->n_mdata / 1e6); + json_float(&json, "mapped_data_fraction", (float)fs->n_mdata / fs->n_data * 100); + json_uint(&json, "properly_paired", fs->n_pair_good); + json_float(&json, "properly_paired_fraction", (float)fs->n_pair_good / fs->n_reads * 100); + json_uint(&json, "read_mate_paired", fs->n_pair_map); + json_float(&json, "read_mate_paired_fraction", (float)fs->n_pair_map / fs->n_reads * 100); + 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", (float)fs->n_dup / fs->n_mapped * 100); + json_int(&json, "mapq_cutoff", f->mapQ_lim); + json_uint(&json, "mapq_reads", fs->n_qual); + json_float(&json, "mapq_reads_fraction_all", (float)fs->n_qual / fs->n_reads * 100); + json_float(&json, "mapq_reads_fraction_mapped", (float)fs->n_qual / fs->n_mapped * 100); + 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", (float)fs->n_tgt / fs->n_reads * 100); + json_float(&json, "target_reads_fraction_mapped", (float)fs->n_tgt / fs->n_mapped * 100); + 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", (float)fs->n_tdata / fs->n_data * 100); + json_float(&json, "target_data_fraction_mapped", (float)fs->n_tdata / fs->n_mdata * 100); + json_uint(&json, "region_length", a->tgt_len); + json_float(&json, "average_depth", (float)fs->n_tdata / a->tgt_len); + json_float(&json, "average_depth_rmdup", (float)fs->n_trmdat / a->tgt_len); + + // 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); + json_uint(&json, "flank_reads", fs->n_flk); + json_float(&json, "flank_reads_fraction_all", (float)fs->n_flk / fs->n_reads * 100); + json_float(&json, "flank_reads_fraction_mapped", (float)fs->n_flk / fs->n_mapped * 100); + json_float(&json, "flank_data_mb", (float)fs->n_fdata / 1e6); + json_float(&json, "flank_data_fraction_all", (float)fs->n_fdata / fs->n_data * 100); + json_float(&json, "flank_data_fraction_mapped", (float)fs->n_fdata / fs->n_mdata * 100); + 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) { int i; @@ -1357,6 +1739,13 @@ 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"); { @@ -1492,6 +1881,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); @@ -1542,6 +1956,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 +1973,7 @@ enum INSERTSIZE, UNCOVER, BAMOUT, + DEPTHRATIO, HELP }; @@ -1568,6 +1986,7 @@ 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}, //{"rmdup", no_argument, NULL, 'd'}, {"help", no_argument, NULL, 'h'}, {"version", no_argument, NULL, 'v'}}; @@ -1632,6 +2051,29 @@ 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; @@ -1730,6 +2172,7 @@ int bamdst(int argc, char *argv[]) if (export_target_bam) bam_close(bamoutfp); freemem(opt.cutoffs); + freemem(opt.ratios); freeall: freemem(export_target_bam); freemem(outdir); From d3ad02be43411b93d7f5d59c0d689bd91ded5aa2 Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 24 Dec 2025 14:02:22 +0800 Subject: [PATCH 02/18] =?UTF-8?q?[update]=20I/O=20=E7=BC=93=E5=86=B2?= =?UTF-8?q?=E4=BC=98=E5=8C=96?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bamdst.c | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/bamdst.c b/bamdst.c index 9450ecc..67250bc 100644 --- a/bamdst.c +++ b/bamdst.c @@ -57,9 +57,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 }; @@ -882,7 +884,11 @@ 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 { @@ -890,6 +896,10 @@ int stat_each_region(loopbams_parameters_t *para, aux_t *a) for (j = 0; j < node->len; ++j) { ksprintf(para->pdepths, "%s\t%d\t0\t0\t0\n", para->name, node->start + j); + } + // 优化：批量写入而不是每行都写 + if (para->pdepths->l > WRITE_BUFFER_SIZE) + { write_buffer_bgzf(para->pdepths, para->fdep); } push_bedreg(para->ucreg, node->start, node->stop); // store uncover region @@ -901,8 +911,11 @@ int stat_each_region(loopbams_parameters_t *para, aux_t *a) 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; } @@ -1519,6 +1532,8 @@ int print_report_json(struct opt_aux *f, aux_t *a, bamflag_t *fs, warnings("Failed to create coverage.report.json"); return -1; } + // 设置更大的文件缓冲区 + setvbuf(fjson, NULL, _IOFBF, WRITE_BUFFER_SIZE); kstring_t json = {0, 0, NULL}; @@ -1701,6 +1716,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(); @@ -1748,6 +1766,8 @@ int print_report(struct opt_aux *f, aux_t *a, bamflag_t *fs) 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 %"); @@ -1808,6 +1828,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); From 110d1aed60d58d0f0c41b70b87abf488d5dac87e Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 24 Dec 2025 14:07:06 +0800 Subject: [PATCH 03/18] =?UTF-8?q?[fix]=20=E4=BF=AE=E5=A4=8D=E9=83=A8?= =?UTF-8?q?=E5=88=86FIXME=E9=97=AE=E9=A2=98?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bamdst.c | 67 ++++++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 58 insertions(+), 9 deletions(-) diff --git a/bamdst.c b/bamdst.c index 67250bc..f85df44 100644 --- a/bamdst.c +++ b/bamdst.c @@ -608,11 +608,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); + int read_ret = bedHand->read(fn, a->h_tgt, 0, 0, &ret); + + // 检查 BED 文件是否被截断或读取失败 + if (read_ret < 0) + { + errabort("Failed to read BED file: %s (file may be truncated or corrupted)", fn); + } + if (zero_based && ret) { warnings("This region is not a standard bed format.\n" @@ -1071,7 +1078,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 @@ -1133,15 +1150,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) { @@ -1434,20 +1456,44 @@ void cntcov_cal_ratio(struct opt_aux *f, struct regcov *cov, count32_t *cnt, uin } } -// FIXME: need improve soon!!! +// 计算深度分布的中位数 +// 使用累积分布查找，时间复杂度 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]; + + if (sum == 0) + return 0; + uint64_t med = sum / 2; uint64_t num = 0; + + // 查找中位数位置 for (i = 0; i < cnt->m; ++i) { num += (uint64_t)cnt->a[i]; if (num >= med) + { + // 对于偶数个元素，返回两个中间值的平均 + if (sum % 2 == 0 && num == med && i + 1 < cnt->m) + { + // 找下一个非零位置 + int j = i + 1; + while (j < cnt->m && cnt->a[j] == 0) + j++; + if (j < cnt->m) + return (float)(i + j) / 2.0f; + } return (float)i; + } } return 0; } @@ -2163,7 +2209,10 @@ int bamdst(int argc, char *argv[]) } aux->ndata = n; } - // FIXME: accpet more than one bam files! + // 多 BAM 文件支持说明： + // 当前实现已支持多个 BAM 文件输入，它们会被顺序处理并合并统计结果 + // 注意：所有 BAM 文件必须使用相同的参考基因组和排序方式 + // 第一个 BAM 文件的 header 会被用于输出 if (export_target_bam) { bamoutfp = bam_open(export_target_bam, "w"); From 32aa999a6b194b442309c74c61e6edac28ac7bed Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 24 Dec 2025 14:13:31 +0800 Subject: [PATCH 04/18] =?UTF-8?q?[fix]=20=E4=BF=AE=E5=A4=8D=E7=BC=96?= =?UTF-8?q?=E8=AF=91=E9=94=99=E8=AF=AF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 2 +- bamdst.c | 15 +++++---------- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/Makefile b/Makefile index 1fb9b73..28edd18 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,5 @@ CC= gcc -CFLAGS= -g -Wall -O2 +CFLAGS= -g -Wall -O2 -std=gnu99 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 diff --git a/bamdst.c b/bamdst.c index f85df44..825321f 100644 --- a/bamdst.c +++ b/bamdst.c @@ -612,12 +612,12 @@ static float coverage_cal(const uint32_t *array, int l) int load_bed_init(char const *fn, aux_t *a) { int ret = 0; - int read_ret = bedHand->read(fn, a->h_tgt, 0, 0, &ret); + bedHand->read(fn, a->h_tgt, 0, 0, &ret); - // 检查 BED 文件是否被截断或读取失败 - if (read_ret < 0) + // 检查 BED 文件是否为空 + if (kh_size(a->h_tgt) == 0) { - errabort("Failed to read BED file: %s (file may be truncated or corrupted)", fn); + errabort("Failed to read BED file: %s (file may be empty or corrupted)", fn); } if (zero_based && ret) @@ -1370,7 +1370,7 @@ uint64_t cntcov_cal_rmdup(struct opt_aux *f, struct regcov *cov, count32_t *cnt, return 0; // 计算各阈值下的覆盖度 - uint64_t below_100 = 0, below_30 = 0, below_10 = 0, below_4 = 0, below_0 = 0; + uint64_t below_100 = 0, below_30 = 0, below_10 = 0, below_4 = 0; for (i = 0; i < cnt->m; ++i) { @@ -1562,11 +1562,6 @@ static void json_float(kstring_t *s, const char *key, float value) ksprintf(s, "\"%s\":%.2f,", key, value); } -static void json_float_array_value(kstring_t *s, float value) -{ - ksprintf(s, "%.2f,", value); -} - /* 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, From dceb23beb3b7391cae8efdd6e1fda0182bee8d43 Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 24 Dec 2025 14:24:10 +0800 Subject: [PATCH 05/18] =?UTF-8?q?[update]=20=E6=A0=BC=E5=BC=8F=E5=8C=96jso?= =?UTF-8?q?n=E8=BE=93=E5=87=BA?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bamdst.c | 33 +++++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/bamdst.c b/bamdst.c index 825321f..2d2da33 100644 --- a/bamdst.c +++ b/bamdst.c @@ -1511,10 +1511,21 @@ float average_cnt(count32_t *cnt) return (float)sum / num; } -/* JSON string builder utilities */ +// 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) @@ -1522,6 +1533,8 @@ 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); } @@ -1539,27 +1552,32 @@ static void json_end_array(kstring_t *s) static void json_key(kstring_t *s, const char *key) { - ksprintf(s, "\"%s\":", key); + json_newline_indent(s); + ksprintf(s, "\"%s\": ", key); } static void json_string(kstring_t *s, const char *key, const char *value) { - ksprintf(s, "\"%s\":\"%s\",", key, value); + json_newline_indent(s); + ksprintf(s, "\"%s\": \"%s\",", key, value); } static void json_int(kstring_t *s, const char *key, int64_t value) { - ksprintf(s, "\"%s\":%" PRId64 ",", key, value); + json_newline_indent(s); + ksprintf(s, "\"%s\": %" PRId64 ",", key, value); } static void json_uint(kstring_t *s, const char *key, uint64_t value) { - ksprintf(s, "\"%s\":%" PRIu64 ",", key, value); + json_newline_indent(s); + ksprintf(s, "\"%s\": %" PRIu64 ",", key, value); } static void json_float(kstring_t *s, const char *key, float value) { - ksprintf(s, "\"%s\":%.2f,", key, value); + json_newline_indent(s); + ksprintf(s, "\"%s\": %.2f,", key, value); } /* Generate JSON format coverage report */ @@ -1576,6 +1594,9 @@ int print_report_json(struct opt_aux *f, aux_t *a, bamflag_t *fs, // 设置更大的文件缓冲区 setvbuf(fjson, NULL, _IOFBF, WRITE_BUFFER_SIZE); + // 重置缩进级别 + json_indent_level = 0; + kstring_t json = {0, 0, NULL}; json_start_object(&json); From c4ad42bc98e0e4ba5af4e9ff210961d047b6519d Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 24 Dec 2025 14:26:20 +0800 Subject: [PATCH 06/18] [release] 1.2.0 --- bamdst.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bamdst.c b/bamdst.c index 2d2da33..81a7a5b 100644 --- a/bamdst.c +++ b/bamdst.c @@ -41,7 +41,7 @@ #include static char const *program_name = "bamdst"; -static char const *Version = "1.1.0"; +static char const *Version = "1.2.0"; /* flank region will be stat in the coverage report file, * this value can be set by -f / --flank */ From 0b9f8090c018bef72409952807470c496f8c5723 Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 13:49:59 +0800 Subject: [PATCH 07/18] =?UTF-8?q?[update]=20=E8=BF=81=E7=A7=BB=E5=88=B0hts?= =?UTF-8?q?lib?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 1. 迁移到htslib 2. 支持cram 3. 帮助信息打印 4. 自动生成输出文件夹 --- Makefile | 33 ++++------- bamdst.c | 163 +++++++++++++++++++++++++++++++++++++------------------ 2 files changed, 122 insertions(+), 74 deletions(-) diff --git a/Makefile b/Makefile index 28edd18..d547fbb 100644 --- a/Makefile +++ b/Makefile @@ -1,10 +1,15 @@ CC= gcc CFLAGS= -g -Wall -O2 -std=gnu99 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 + +# htslib support - use pkg-config if available, otherwise use default paths +HTSLIB_CFLAGS := $(shell pkg-config --cflags htslib 2>/dev/null || echo "") +HTSLIB_LIBS := $(shell pkg-config --libs htslib 2>/dev/null || echo "-lhts") + +# Objects for output file writing (bgzf for depth.tsv.gz etc.) +LOBJS= bgzf.o kstring.o bedutil.o commons.o PROG= bamdst -INCLUDES= -Isamlib/ -I. -SUBDIRS= . samlib +INCLUDES= -I. $(HTSLIB_CFLAGS) LIBPATH= -L. .SUFFIXES:.c .o @@ -15,15 +20,15 @@ LIBPATH= -L. all:clean $(PROG) -.PHONY:all clean +.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 +bamdst:lib + $(CC) $(CFLAGS) -o $@ $(LDFLAGS) bamdst.c $(LIBPATH) $(INCLUDES) -lm -lbam $(HTSLIB_LIBS) -lz -lpthread bgzf.o:bgzf.c bgzf.h $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) bgzf.c -o $@ @@ -31,21 +36,6 @@ bgzf.o:bgzf.c bgzf.h 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 $@ @@ -54,3 +44,4 @@ bedutil.o:bedutil.c bedutil.h clean: rm -fr gmon.out *.o a.out *.exe *.dSYM $(PROG) *~ *.a target.dep *.plot *.report *.tsv.gz uncover.bed + diff --git a/bamdst.c b/bamdst.c index 81a7a5b..e34f229 100644 --- a/bamdst.c +++ b/bamdst.c @@ -1,5 +1,6 @@ /* The MIT License + Copyright (c) 2026 pzweuj - Ported to HTSlib Copyright (c) 2022, 2023, 2024 Authors Copyright (c) 2013-2014 Beijing Genomics Institution (BGI) @@ -29,12 +30,15 @@ #include "commons.h" #include "count.h" -// bam.h and sam_header.h are standard header from samtools -#include "bam.h" -#include "sam_header.h" +// htslib for BAM/CRAM/SAM file reading +#include +#include +#include + +// bgzf for writing tabix-able depth.gz file (keep local implementation) +#include "bgzf.h" // 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" @@ -105,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() @@ -161,6 +163,8 @@ struct opt_aux int num_ratios; // 比例数量 int max_ratios; // 最大比例数量 float *ratios; // 比例数组 (如 0.1, 0.2, 0.5) + // 新增：参考基因组路径（用于 CRAM 文件） + char *reference; // 参考基因组 FASTA 文件路径 }; // 一个函数来初始化 opt_aux 结构体 @@ -197,6 +201,9 @@ struct opt_aux init_opt_aux() opt.ratios[0] = 0.2f; opt.ratios[1] = 0.5f; + // 初始化参考基因组路径 + opt.reference = NULL; + return opt; } @@ -309,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 */ @@ -334,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); @@ -357,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); */ @@ -465,17 +472,19 @@ 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\ @@ -489,6 +498,10 @@ Optional parameters:\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\ @@ -503,9 +516,11 @@ Optional parameters:\n\ - uncover.bed the bad covered or uncovered region in the probe file\n\ \n\ * New features in this version:\n\ + - CRAM format support: use -T to specify reference genome\n\ - Rmdup coverage statistics: coverage based on deduplicated depth\n\ - Ratio-based coverage: coverage at ratios of average depth (e.g., 0.2x, 0.5x)\n\ - JSON output: coverage.report.json for programmatic access\n\ + - Auto directory creation: output directory created if not exists\n\ \n\ * About depth.tsv.gz:\n\ * There are five columns in this file, including chromosome, position, raw\n\ @@ -999,7 +1014,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"); @@ -1008,25 +1023,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 @@ -1115,7 +1130,7 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) } 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 @@ -1190,7 +1205,7 @@ int load_bamfiles(struct opt_aux *f, aux_t *a, bamflag_t *fs) // endcore: } bam_destroy1(b); - bgzf_close(a->data[i]); + hts_close(a->data[i]); } while (para->tgt_node) { @@ -2058,6 +2073,7 @@ enum UNCOVER, BAMOUT, DEPTHRATIO, + REFERENCE, HELP }; @@ -2071,6 +2087,7 @@ static struct option const long_opts[] = {{"outdir", required_argument, NULL, 'o {"uncover", required_argument, NULL, UNCOVER}, {"bamout", required_argument, NULL, BAMOUT}, {"depthratio", required_argument, NULL, DEPTHRATIO}, + {"reference", required_argument, NULL, 'T'}, //{"rmdup", no_argument, NULL, 'd'}, {"help", no_argument, NULL, 'h'}, {"version", no_argument, NULL, 'v'}}; @@ -2087,7 +2104,7 @@ int bamdst(int argc, char *argv[]) // 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) { @@ -2161,6 +2178,10 @@ int bamdst(int argc, char *argv[]) case 'q': opt.mapQ_lim = atoi(optarg); break; + case 'T': + // 参考基因组路径（用于 CRAM 文件） + opt.reference = strdup(optarg); + break; case 'h': usage(1); break; @@ -2176,7 +2197,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); @@ -2185,56 +2206,91 @@ 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"); + // 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]); + + // 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; } - // 多 BAM 文件支持说明： - // 当前实现已支持多个 BAM 文件输入，它们会被顺序处理并合并统计结果 - // 注意：所有 BAM 文件必须使用相同的参考基因组和排序方式 - // 第一个 BAM 文件的 header 会被用于输出 + // 多 BAM/CRAM 文件支持说明： + // 当前实现已支持多个 BAM/CRAM 文件输入，它们会被顺序处理并合并统计结果 + // 注意：所有文件必须使用相同的参考基因组和排序方式 + // 第一个文件的 header 会被用于输出 if (export_target_bam) { - bamoutfp = bam_open(export_target_bam, "w"); + bamoutfp = bgzf_open(export_target_bam, "w"); if (bamoutfp == NULL) errabort("%s : %s", export_target_bam, strerror(errno)); - bam_header_write(bamoutfp, aux->h); + // Write BAM header using htslib compatible format + if (bam_hdr_write(bamoutfp, aux->h) < 0) + errabort("Failed to write header to %s", export_target_bam); } h_chrlength_init(); header2chrhash(aux->h); @@ -2248,7 +2304,7 @@ int bamdst(int argc, char *argv[]) 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); @@ -2257,9 +2313,10 @@ int bamdst(int argc, char *argv[]) freemem(opt.inputs[i]); freemem(opt.inputs); if (export_target_bam) - bam_close(bamoutfp); + bgzf_close(bamoutfp); freemem(opt.cutoffs); freemem(opt.ratios); + freemem(opt.reference); freeall: freemem(export_target_bam); freemem(outdir); From d7036676fd862c2cac0243082b0d0c01ae8e734b Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 14:04:03 +0800 Subject: [PATCH 08/18] [release] 2.0.0 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 迁移到HTSLib，增加多线程支持和CRAM支持 --- Makefile | 2 +- README.md | 35 +- bamdst.c | 55 ++- samlib/bam.c | 474 -------------------------- samlib/bam.h | 795 ------------------------------------------- samlib/bam_aux.c | 217 ------------ samlib/bam_endian.h | 42 --- samlib/bam_import.c | 489 -------------------------- samlib/bam_index.c | 726 --------------------------------------- samlib/sam_header.c | 810 -------------------------------------------- samlib/sam_header.h | 48 --- 11 files changed, 76 insertions(+), 3617 deletions(-) delete mode 100644 samlib/bam.c delete mode 100644 samlib/bam.h delete mode 100644 samlib/bam_aux.c delete mode 100644 samlib/bam_endian.h delete mode 100644 samlib/bam_import.c delete mode 100644 samlib/bam_index.c delete mode 100644 samlib/sam_header.c delete mode 100644 samlib/sam_header.h diff --git a/Makefile b/Makefile index d547fbb..3628b92 100644 --- a/Makefile +++ b/Makefile @@ -28,7 +28,7 @@ libbam.a:$(LOBJS) $(AR) -csru $@ $(LOBJS) bamdst:lib - $(CC) $(CFLAGS) -o $@ $(LDFLAGS) bamdst.c $(LIBPATH) $(INCLUDES) -lm -lbam $(HTSLIB_LIBS) -lz -lpthread + $(CC) $(CFLAGS) -o $@ $(LDFLAGS) bamdst.c $(LIBPATH) $(INCLUDES) -lm $(HTSLIB_LIBS) -lz -lpthread bgzf.o:bgzf.c bgzf.h $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) bgzf.c -o $@ diff --git a/README.md b/README.md index 0930894..60f405f 100644 --- a/README.md +++ b/README.md @@ -1,30 +1,44 @@ # bamdst -- a BAM Depth Stat Tool -Bamdst is a lightweight tool to stat the depth coverage of target regions of bam file(s). +Bamdst is a lightweight tool to stat the depth coverage of target regions of bam file(s). Bam file(s) should be properly sorted, and the probe file (bed file) and the output dir must be specified in the first place. -## Install +## Dependencies + +- **htslib** - Required for BAM/CRAM/SAM file reading and writing +- **zlib** - Required for compression +- **pthread** - Required for multi-threading support + +On Debian/Ubuntu: +```bash +apt-get install libhts-dev zlib1g-dev ``` -git clone https://github.com/shiquan/bamdst -cd bamdst -make + +On macOS with Homebrew: +```bash +brew install htslib ``` -Or by conda. Thanks @xdgene +## Install ``` -conda install xdgene::bamdst +git clone https://github.com/pzweuj/bamdst +cd bamdst +make ``` - ## USAGE Normal: bamdst -p -o ./ in1.bam +With multi-threading (8 threads): + + bamdst -p -o ./ --threads 8 in1.bam + Pipeline mode: samtools view in1.bam -u | bamdst -p x.bed -o ./ - @@ -79,6 +93,11 @@ Use comma-separated values like "0.1,0.2,0.5". Default is "0.2,0.5". This calculates coverage at positions with depth > (ratio × average_depth). +--threads [num] + +set the number of threads for BAM/CRAM I/O operations. Uses htslib's native +multi-threading for faster decompression and compression. Default is 0 (single-threaded mode). + ## OUTPUT FILES Eight files will be created in the output direction. There are: diff --git a/bamdst.c b/bamdst.c index e34f229..57cce1d 100644 --- a/bamdst.c +++ b/bamdst.c @@ -45,7 +45,7 @@ #include static char const *program_name = "bamdst"; -static char const *Version = "1.2.0"; +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 */ @@ -77,7 +77,7 @@ static const int WRITE_BUFFER_SIZE = 65536; /* 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) { @@ -165,6 +165,8 @@ struct opt_aux float *ratios; // 比例数组 (如 0.1, 0.2, 0.5) // 新增：参考基因组路径（用于 CRAM 文件） char *reference; // 参考基因组 FASTA 文件路径 + // 新增：线程数 + int nthreads; // htslib 多线程压缩/解压缩线程数 }; // 一个函数来初始化 opt_aux 结构体 @@ -203,7 +205,10 @@ struct opt_aux init_opt_aux() // 初始化参考基因组路径 opt.reference = NULL; - + + // 初始化线程数（0 = 单线程模式，不使用多线程） + opt.nthreads = 0; + return opt; } @@ -493,6 +498,7 @@ Optional parameters:\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"); @@ -2074,6 +2080,7 @@ enum BAMOUT, DEPTHRATIO, REFERENCE, + THREADS, HELP }; @@ -2088,6 +2095,7 @@ static struct option const long_opts[] = {{"outdir", required_argument, NULL, 'o {"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'}}; @@ -2182,6 +2190,14 @@ int bamdst(int argc, char *argv[]) // 参考基因组路径（用于 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; @@ -2215,6 +2231,14 @@ int bamdst(int argc, char *argv[]) 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) { @@ -2248,7 +2272,16 @@ int bamdst(int argc, char *argv[]) } if (aux->data[i] == NULL) errabort("%s: %s", argv[optind + i], strerror(errno)); - + + // 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) @@ -2285,11 +2318,19 @@ int bamdst(int argc, char *argv[]) // 第一个文件的 header 会被用于输出 if (export_target_bam) { - bamoutfp = bgzf_open(export_target_bam, "w"); + bamoutfp = hts_open(export_target_bam, "wb"); if (bamoutfp == NULL) errabort("%s : %s", export_target_bam, strerror(errno)); + // 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 (bam_hdr_write(bamoutfp, aux->h) < 0) + if (sam_hdr_write(bamoutfp, aux->h) < 0) errabort("Failed to write header to %s", export_target_bam); } h_chrlength_init(); @@ -2313,7 +2354,7 @@ int bamdst(int argc, char *argv[]) freemem(opt.inputs[i]); freemem(opt.inputs); if (export_target_bam) - bgzf_close(bamoutfp); + hts_close(bamoutfp); freemem(opt.cutoffs); freemem(opt.ratios); freemem(opt.reference); 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 Date: Tue, 6 Jan 2026 14:13:36 +0800 Subject: [PATCH 09/18] =?UTF-8?q?[fix]=20=E8=B0=83=E7=94=A8=E9=94=99?= =?UTF-8?q?=E8=AF=AF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bamdst.c | 6 +++--- bgzf.h | 5 +++++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/bamdst.c b/bamdst.c index 57cce1d..2dc7d63 100644 --- a/bamdst.c +++ b/bamdst.c @@ -756,8 +756,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_cigar2pos(cigar); if (end >= tmp->start) { @@ -1204,7 +1204,7 @@ 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, aux->h, b); fs->n_tgt++; } diff --git a/bgzf.h b/bgzf.h index 29fe0e5..6e88a8a 100644 --- a/bgzf.h +++ b/bgzf.h @@ -32,6 +32,10 @@ #include #include +/* Use htslib's BGZF if available, otherwise define our own */ +#ifdef HTS_READAHEAD_SIZE +#include +#else #define BGZF_BLOCK_SIZE 0x10000 #define BGZF_MAX_BLOCK_SIZE 0x10000 @@ -49,6 +53,7 @@ typedef struct { void *cache; // a pointer to a hash table void *fp; // actual file handler; FILE* on writing; FILE* or knetFile* on reading } BGZF; +#endif #ifndef KSTRING_T #define KSTRING_T kstring_t From 5030967d5e784ef50c19481fe6fe99c10bf0f0a7 Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 14:15:52 +0800 Subject: [PATCH 10/18] =?UTF-8?q?[fix]=20=E8=B0=83=E7=94=A8=E4=BF=AE?= =?UTF-8?q?=E5=A4=8D?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bamdst.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/bamdst.c b/bamdst.c index 2dc7d63..dea0ff8 100644 --- a/bamdst.c +++ b/bamdst.c @@ -30,14 +30,14 @@ #include "commons.h" #include "count.h" +// bgzf for writing tabix-able depth.gz file (must be before htslib to avoid BGZF conflict) +#include "bgzf.h" + // htslib for BAM/CRAM/SAM file reading #include #include #include -// bgzf for writing tabix-able depth.gz file (keep local implementation) -#include "bgzf.h" - // khash, kstring and knetfile are standard utils of klib #include "khash.h" #include "knetfile.h" @@ -757,7 +757,7 @@ int readcore(struct depnode *header, bam1_t const *b, cntstat_t state) if (isNull(tmp)) return 0; uint32_t *cigar = bam_get_cigar(b); - uint32_t end = bam_cigar2pos(cigar); + uint32_t end = bam_endpos(b); if (end >= tmp->start) { @@ -1204,7 +1204,7 @@ 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) - sam_write1(bamoutfp, aux->h, b); + sam_write1(bamoutfp, a->h, b); fs->n_tgt++; } From 4ac3c4a5058a948df0527e7ca189cfc7730e41bc Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 14:18:46 +0800 Subject: [PATCH 11/18] =?UTF-8?q?[fix]=20=E5=BC=83=E7=94=A8=E6=9C=AC?= =?UTF-8?q?=E5=9C=B0bgzf?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 9 +- bamdst.c | 4 +- bgzf.c | 555 ------------------------------------------------------- bgzf.h | 201 -------------------- 4 files changed, 5 insertions(+), 764 deletions(-) delete mode 100644 bgzf.c delete mode 100644 bgzf.h diff --git a/Makefile b/Makefile index 3628b92..6ab9fb6 100644 --- a/Makefile +++ b/Makefile @@ -1,13 +1,13 @@ CC= gcc CFLAGS= -g -Wall -O2 -std=gnu99 -DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -DBGZF_CACHE +DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE # htslib support - use pkg-config if available, otherwise use default paths HTSLIB_CFLAGS := $(shell pkg-config --cflags htslib 2>/dev/null || echo "") HTSLIB_LIBS := $(shell pkg-config --libs htslib 2>/dev/null || echo "-lhts") -# Objects for output file writing (bgzf for depth.tsv.gz etc.) -LOBJS= bgzf.o kstring.o bedutil.o commons.o +# Objects for output file writing (bgzf is now from htslib) +LOBJS= kstring.o bedutil.o commons.o PROG= bamdst INCLUDES= -I. $(HTSLIB_CFLAGS) LIBPATH= -L. @@ -30,9 +30,6 @@ libbam.a:$(LOBJS) bamdst:lib $(CC) $(CFLAGS) -o $@ $(LDFLAGS) bamdst.c $(LIBPATH) $(INCLUDES) -lm $(HTSLIB_LIBS) -lz -lpthread -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 $@ diff --git a/bamdst.c b/bamdst.c index dea0ff8..c2457d3 100644 --- a/bamdst.c +++ b/bamdst.c @@ -30,8 +30,8 @@ #include "commons.h" #include "count.h" -// bgzf for writing tabix-able depth.gz file (must be before htslib to avoid BGZF conflict) -#include "bgzf.h" +// bgzf for writing tabix-able depth.gz file +#include // htslib for BAM/CRAM/SAM file reading #include 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 6e88a8a..0000000 --- a/bgzf.h +++ /dev/null @@ -1,201 +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 - -/* Use htslib's BGZF if available, otherwise define our own */ -#ifdef HTS_READAHEAD_SIZE -#include -#else -#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; -#endif - -#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 From 6566a7adc1c52d7c8ac7b23fe2e600ce8040088b Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 14:21:39 +0800 Subject: [PATCH 12/18] =?UTF-8?q?[fix]=20=E4=BF=AE=E6=AD=A3makefile?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 25 +++++-------------------- 1 file changed, 5 insertions(+), 20 deletions(-) diff --git a/Makefile b/Makefile index 6ab9fb6..d6c027c 100644 --- a/Makefile +++ b/Makefile @@ -6,11 +6,10 @@ DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE HTSLIB_CFLAGS := $(shell pkg-config --cflags htslib 2>/dev/null || echo "") HTSLIB_LIBS := $(shell pkg-config --libs htslib 2>/dev/null || echo "-lhts") -# Objects for output file writing (bgzf is now from htslib) -LOBJS= kstring.o bedutil.o commons.o +# Source files +SOURCES= bamdst.c kstring.c bedutil.c commons.c PROG= bamdst INCLUDES= -I. $(HTSLIB_CFLAGS) -LIBPATH= -L. .SUFFIXES:.c .o .PHONY: all @@ -18,26 +17,12 @@ LIBPATH= -L. .c.o: $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) $< -o $@ -all:clean $(PROG) +all:clean $(PROG) .PHONY:all clean -lib:libbam.a - -libbam.a:$(LOBJS) - $(AR) -csru $@ $(LOBJS) - -bamdst:lib - $(CC) $(CFLAGS) -o $@ $(LDFLAGS) bamdst.c $(LIBPATH) $(INCLUDES) -lm $(HTSLIB_LIBS) -lz -lpthread - -kstring.o:kstring.c kstring.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) kstring.c -o $@ - -commons.o:commons.c commons.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) commons.c -o $@ - -bedutil.o:bedutil.c bedutil.h - $(CC) -c $(CFLAGS) $(DFLAGS) $(INCLUDES) bedutil.c -o $@ +bamdst: $(SOURCES:.c=.o) + $(CC) $(CFLAGS) -o $@ $(SOURCES:.c=.o) $(INCLUDES) -lm $(HTSLIB_LIBS) -lz -lpthread clean: rm -fr gmon.out *.o a.out *.exe *.dSYM $(PROG) *~ *.a target.dep *.plot *.report *.tsv.gz uncover.bed From 63c4de20fdb7f33c0335d4904b5a1a39b0359656 Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 14:31:43 +0800 Subject: [PATCH 13/18] =?UTF-8?q?[=E5=88=86=E6=94=AF]=20=E4=BB=8Ebamdst?= =?UTF-8?q?=E5=88=86=E6=94=AF=E5=88=B0xamdst?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 4 ++-- README.md | 18 ++++++++++-------- bamdst.c => xamdst.c | 15 ++++----------- 3 files changed, 16 insertions(+), 21 deletions(-) rename bamdst.c => xamdst.c (99%) diff --git a/Makefile b/Makefile index d6c027c..899d865 100644 --- a/Makefile +++ b/Makefile @@ -7,8 +7,8 @@ HTSLIB_CFLAGS := $(shell pkg-config --cflags htslib 2>/dev/null || echo "") HTSLIB_LIBS := $(shell pkg-config --libs htslib 2>/dev/null || echo "-lhts") # Source files -SOURCES= bamdst.c kstring.c bedutil.c commons.c -PROG= bamdst +SOURCES= xamdst.c kstring.c bedutil.c commons.c +PROG= xamdst INCLUDES= -I. $(HTSLIB_CFLAGS) .SUFFIXES:.c .o diff --git a/README.md b/README.md index 60f405f..8efee80 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,10 @@ -# bamdst -- a BAM Depth Stat Tool +# xamdst -- a BAM/CRAM Depth Stat Tool -Bamdst is a lightweight tool to stat the depth coverage of target regions of bam file(s). +This project is a hard fork of bamdst by pzweuj. It has been refactored to use HTSlib instead of the legacy samlib/bam.h, bringing support for multi-threading and CRAM format. -Bam file(s) should be properly sorted, and the probe file (bed file) and the output dir +xamdst is a lightweight tool to stat the depth coverage of target regions of BAM/CRAM/SAM file(s). + +Input files should be properly sorted, and the probe file (bed file) and the output dir must be specified in the first place. @@ -24,8 +26,8 @@ brew install htslib ## Install ``` -git clone https://github.com/pzweuj/bamdst -cd bamdst +git clone https://github.com/pzweuj/xamdst +cd xamdst make ``` @@ -33,15 +35,15 @@ make Normal: - bamdst -p -o ./ in1.bam + xamdst -p -o ./ in1.bam With multi-threading (8 threads): - bamdst -p -o ./ --threads 8 in1.bam + xamdst -p -o ./ --threads 8 in1.bam Pipeline mode: - samtools view in1.bam -u | bamdst -p x.bed -o ./ - + samtools view in1.bam -u | xamdst -p x.bed -o ./ - ## PARAMETERS diff --git a/bamdst.c b/xamdst.c similarity index 99% rename from bamdst.c rename to xamdst.c index c2457d3..8e99e08 100644 --- a/bamdst.c +++ b/xamdst.c @@ -44,7 +44,7 @@ #include "kstring.h" #include -static char const *program_name = "bamdst"; +static char const *program_name = "xamdst"; static char const *Version = "2.0.0"; /* flank region will be stat in the coverage report file, @@ -521,13 +521,6 @@ Optional parameters:\n\ - depth.tsv.gz raw depth, rmdup depth, coverage depth of each position\n\ - uncover.bed the bad covered or uncovered region in the probe file\n\ \n\ -* New features in this version:\n\ - - CRAM format support: use -T to specify reference genome\n\ - - Rmdup coverage statistics: coverage based on deduplicated depth\n\ - - Ratio-based coverage: coverage at ratios of average depth (e.g., 0.2x, 0.5x)\n\ - - JSON output: coverage.report.json for programmatic access\n\ - - Auto directory creation: output directory created if not exists\n\ -\n\ * About depth.tsv.gz:\n\ * There are five columns in this file, including chromosome, position, raw\n\ * depth, rmdep depth, coverage depth\n\ @@ -541,7 +534,7 @@ Optional parameters:\n\ "); puts("============\n"); puts(" HOMEPAGE: \n\ - https://github.com/shiquan/bamdst\n"); + https://github.com/pzweuj/xamdst\n"); } exit(EXIT_SUCCESS); } @@ -2105,7 +2098,7 @@ 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; @@ -2367,5 +2360,5 @@ int bamdst(int argc, char *argv[]) /* main */ int main(int argc, char *argv[]) { - return bamdst(argc, argv); + return xamdst(argc, argv); } From 34be3c9cfcc0672c41ed6a59f00ed75173aa8e9f Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 14:33:28 +0800 Subject: [PATCH 14/18] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 8efee80..dbf8cae 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # xamdst -- a BAM/CRAM Depth Stat Tool -This project is a hard fork of bamdst by pzweuj. It has been refactored to use HTSlib instead of the legacy samlib/bam.h, bringing support for multi-threading and CRAM format. +This project is a hard fork of [bamdst by shiquan](https://github.com/shiquan/bamdst). It has been refactored to use HTSlib instead of the legacy samlib/bam.h, bringing support for multi-threading and CRAM format. xamdst is a lightweight tool to stat the depth coverage of target regions of BAM/CRAM/SAM file(s). From 6b44c0f7a7f104b31501dd5c75ec3b7d3aabc382 Mon Sep 17 00:00:00 2001 From: EuJ Date: Tue, 6 Jan 2026 15:57:23 +0800 Subject: [PATCH 15/18] Update Makefile --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 899d865..e4905be 100644 --- a/Makefile +++ b/Makefile @@ -21,7 +21,7 @@ all:clean $(PROG) .PHONY:all clean -bamdst: $(SOURCES:.c=.o) +xamdst: $(SOURCES:.c=.o) $(CC) $(CFLAGS) -o $@ $(SOURCES:.c=.o) $(INCLUDES) -lm $(HTSLIB_LIBS) -lz -lpthread clean: From f11fa7f4b198972cf8bcb3a0f848ceaf1cceba50 Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 7 Jan 2026 10:05:28 +0800 Subject: [PATCH 16/18] Include Docker pull command in README Add Docker instructions for pulling the mapping image --- README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/README.md b/README.md index dbf8cae..3f4133d 100644 --- a/README.md +++ b/README.md @@ -31,6 +31,11 @@ cd xamdst make ``` +or Docker +``` +docker pull ghcr.io/pzweuj/mapping:2026Jan +``` + ## USAGE Normal: From 55b34852bc156fbeb15d4138eb1e9af9d19e3e91 Mon Sep 17 00:00:00 2001 From: EuJ Date: Wed, 8 Jul 2026 09:59:19 +0800 Subject: [PATCH 17/18] =?UTF-8?q?=E6=95=88=E7=8E=87=E4=BC=98=E5=8C=96?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 67 ++++++++++++++--- count.h | 22 ++++-- xamdst.c | 217 ++++++++++++++++++++++++++++++++++--------------------- 3 files changed, 208 insertions(+), 98 deletions(-) diff --git a/Makefile b/Makefile index e4905be..4e35554 100644 --- a/Makefile +++ b/Makefile @@ -1,29 +1,74 @@ CC= gcc -CFLAGS= -g -Wall -O2 -std=gnu99 +CFLAGS= -g -Wall -O2 -std=gnu99 DFLAGS= -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -# htslib support - use pkg-config if available, otherwise use default paths -HTSLIB_CFLAGS := $(shell pkg-config --cflags htslib 2>/dev/null || echo "") -HTSLIB_LIBS := $(shell pkg-config --libs htslib 2>/dev/null || echo "-lhts") - # 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 - -xamdst: $(SOURCES:.c=.o) - $(CC) $(CFLAGS) -o $@ $(SOURCES:.c=.o) $(INCLUDES) -lm $(HTSLIB_LIBS) -lz -lpthread +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/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/xamdst.c b/xamdst.c index 8e99e08..74de425 100644 --- a/xamdst.c +++ b/xamdst.c @@ -218,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; @@ -240,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) { @@ -252,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; @@ -271,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); \ @@ -295,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; @@ -592,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; } @@ -696,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)) @@ -728,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 @@ -759,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) { @@ -774,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++; } } @@ -856,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 @@ -874,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); @@ -913,10 +938,15 @@ int stat_each_region(loopbams_parameters_t *para, aux_t *a) } 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) { - ksprintf(para->pdepths, "%s\t%d\t0\t0\t0\n", para->name, node->start + j); + kput_depth_line(para->pdepths, para->name, name_len, (int)(node->start + j), 0, 0, 0); } // 优化：批量写入而不是每行都写 if (para->pdepths->l > WRITE_BUFFER_SIZE) @@ -924,8 +954,9 @@ int stat_each_region(loopbams_parameters_t *para, aux_t *a) 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); @@ -968,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 } @@ -1123,7 +1157,8 @@ 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); } } @@ -1213,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); @@ -1311,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) @@ -1487,29 +1525,31 @@ float median_cnt(count32_t *cnt) if (sum == 0) return 0; - uint64_t med = sum / 2; + // 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) + if (v1 < 0 && num >= lower) + v1 = i; + if (v2 < 0 && num >= upper) { - // 对于偶数个元素，返回两个中间值的平均 - if (sum % 2 == 0 && num == med && i + 1 < cnt->m) - { - // 找下一个非零位置 - int j = i + 1; - while (j < cnt->m && cnt->a[j] == 0) - j++; - if (j < cnt->m) - return (float)(i + j) / 2.0f; - } - return (float)i; + 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) @@ -1522,7 +1562,7 @@ 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 格式化输出 - 缩进级别跟踪 @@ -1594,6 +1634,13 @@ static void json_float(kstring_t *s, const char *key, float value) 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, @@ -1633,13 +1680,13 @@ int print_report_json(struct opt_aux *f, aux_t *a, bamflag_t *fs, 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", (float)fs->n_mapped / fs->n_reads * 100); + 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", (float)fs->n_mdata / fs->n_data * 100); + 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", (float)fs->n_pair_good / fs->n_reads * 100); + 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", (float)fs->n_pair_map / fs->n_reads * 100); + 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); @@ -1649,11 +1696,11 @@ int print_report_json(struct opt_aux *f, aux_t *a, bamflag_t *fs, 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", (float)fs->n_dup / fs->n_mapped * 100); + 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", (float)fs->n_qual / fs->n_reads * 100); - json_float(&json, "mapq_reads_fraction_mapped", (float)fs->n_qual / fs->n_mapped * 100); + 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); @@ -1669,15 +1716,15 @@ int print_report_json(struct opt_aux *f, aux_t *a, bamflag_t *fs, json_key(&json, "target"); json_start_object(&json); json_uint(&json, "target_reads", fs->n_tgt); - json_float(&json, "target_reads_fraction_all", (float)fs->n_tgt / fs->n_reads * 100); - json_float(&json, "target_reads_fraction_mapped", (float)fs->n_tgt / fs->n_mapped * 100); + 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", (float)fs->n_tdata / fs->n_data * 100); - json_float(&json, "target_data_fraction_mapped", (float)fs->n_tdata / fs->n_mdata * 100); + 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); - json_float(&json, "average_depth_rmdup", (float)fs->n_trmdat / 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"); @@ -1756,13 +1803,13 @@ int print_report_json(struct opt_aux *f, aux_t *a, bamflag_t *fs, 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); + 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", (float)fs->n_flk / fs->n_reads * 100); - json_float(&json, "flank_reads_fraction_mapped", (float)fs->n_flk / fs->n_mapped * 100); + 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", (float)fs->n_fdata / fs->n_data * 100); - json_float(&json, "flank_data_fraction_mapped", (float)fs->n_fdata / fs->n_mdata * 100); + 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); @@ -1919,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); @@ -1938,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); @@ -2026,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); @@ -2331,13 +2378,19 @@ int xamdst(int argc, char *argv[]) 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 = sam_hdr_nref(aux->h); struct bamflag fs = {}; load_bamfiles(&opt, aux, &fs); From f6f99ad2bb1b2a3c7c0591dbe8e46296be67ce23 Mon Sep 17 00:00:00 2001 From: EuJ Date: Thu, 9 Jul 2026 10:23:41 +0800 Subject: [PATCH 18/18] Update README.md --- README.md | 314 +++++++++++++++--------------------------------------- 1 file changed, 88 insertions(+), 226 deletions(-) diff --git a/README.md b/README.md index 3f4133d..71cb5f2 100644 --- a/README.md +++ b/README.md @@ -1,253 +1,115 @@ -# xamdst -- a BAM/CRAM Depth Stat Tool +# xamdst -This project is a hard fork of [bamdst by shiquan](https://github.com/shiquan/bamdst). It has been refactored to use HTSlib instead of the legacy samlib/bam.h, bringing support for multi-threading and CRAM format. +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. -xamdst is a lightweight tool to stat the depth coverage of target regions of BAM/CRAM/SAM file(s). +Input files must be **coordinate-sorted**. The probe file and output directory are +mandatory. -Input files should be properly sorted, and the probe file (bed file) and the output dir +## Features -must be specified in the first place. +- 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 ## Dependencies -- **htslib** - Required for BAM/CRAM/SAM file reading and writing -- **zlib** - Required for compression -- **pthread** - Required for multi-threading support +- **htslib**, **zlib**, **pthread** -On Debian/Ubuntu: ```bash +# Debian/Ubuntu apt-get install libhts-dev zlib1g-dev -``` - -On macOS with Homebrew: -```bash +# macOS (Homebrew) brew install htslib ``` ## Install -``` + +```bash git clone https://github.com/pzweuj/xamdst cd xamdst make ``` -or Docker -``` +Or use the Docker image: + +```bash docker pull ghcr.io/pzweuj/mapping:2026Jan ``` -## USAGE - -Normal: - - xamdst -p -o ./ in1.bam - -With multi-threading (8 threads): - - xamdst -p -o ./ --threads 8 in1.bam - -Pipeline mode: - - samtools view in1.bam -u | xamdst -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. - ---depthratio [ratios] - -specify depth ratios for coverage calculation relative to average depth. - -Use comma-separated values like "0.1,0.2,0.5". Default is "0.2,0.5". - -This calculates coverage at positions with depth > (ratio × average_depth). - ---threads [num] - -set the number of threads for BAM/CRAM I/O operations. Uses htslib's native -multi-threading for faster decompression and compression. Default is 0 (single-threaded mode). +## Usage -## OUTPUT FILES - -Eight files will be created in the output direction. There are: - --**coverage.report** - --**coverage.report.json** (NEW: JSON format for programmatic access) - --**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 - -the raw depth with consider of deletion region, so this value - -should be equal to or greated than the raw depth. We usw raw depth - -to stat the coverage information in the coverage.report file. If - -you want use rmdup depth to calculate the coverage, please use - -the parameter "--use_rmdup". - -**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: