bam_plp_autoの使用例

bam_plp_auto.c

#include <stdio.h>
#include <stdlib.h>
#include <htslib/sam.h>

typedef struct {
    samFile *fp;
    bam_hdr_t *h;
} bam_file_handle;

static int process_reads(void *data, bam1_t *b) {
    bam_file_handle *aux = (bam_file_handle *)data;
    int ret;
    ret = sam_read1(aux->fp, aux->h, b);

    return ret;
}

int main(int argc, char *argv[]) {
    if (argc < 2) {
        fprintf(stderr, "Usage: %s <bam_file> <min_mapQ>\n", argv[0]);
        return 1;
    }

    samFile *fp = sam_open(argv[1], "r");
    if (fp == NULL) {
        fprintf(stderr, "Error opening BAM file.\n");
        return 1;
    }

    bam_hdr_t *hdr = sam_hdr_read(fp);
    if (hdr == NULL) {
        fprintf(stderr, "Error reading header from BAM file.\n");
        sam_close(fp);
        return 1;
    }

    bam_file_handle aux = {fp, hdr};

    bam_plp_t plp = bam_plp_init(process_reads, &aux);
    if (plp == NULL) {
        fprintf(stderr, "Failed to initialize pileup.\n");
        bam_hdr_destroy(hdr);
        sam_close(fp);
        return 1;
    }

    const bam_pileup1_t *p;
    int tid, pos, n_plp;
    while ((p = bam_plp_auto(plp, &tid, &pos, &n_plp)) != NULL) {
        printf("Position %d on %s: %d reads", pos, hdr->target_name[tid], n_plp);
        for(int i = 0; i < n_plp; i++) {
            const bam_pileup1_t *pp = p + i;
            // puts read names
            printf(" %s\t", bam_get_qname(pp->b));
            printf(" qpos=%d, indel=%d, level=%d, is_del=%d, is_head=%d, is_tail=%d, is_refskip=%d, aux=%u\n",
                     pp->qpos, pp->indel, pp->level, pp->is_del, pp->is_head, pp->is_tail, pp->is_refskip, pp->aux);
        }
    }

    // リソースの解放
    bam_plp_destroy(plp);
    bam_hdr_destroy(hdr);
    sam_close(fp);

    return 0;
}