pileup sample by ChatGPT

t.c

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

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

    // BAMファイルを開く
    samFile *fp = sam_open(argv[1], "r");
    if (fp == NULL) {
        fprintf(stderr, "Failed to open BAM file: %s\n", argv[1]);
        return 1;
    }

    // BAMヘッダーを読み込む
    bam_hdr_t *hdr = sam_hdr_read(fp);
    if (hdr == NULL) {
        fprintf(stderr, "Failed to read BAM header.\n");
        sam_close(fp);
        return 1;
    }

    // Pileupイテレータを初期化
    bam_plp_t plp = bam_plp_init(NULL, NULL);  // コールバックなしで初期化
    if (plp == NULL) {
        fprintf(stderr, "Failed to initialize pileup.\n");
        bam_hdr_destroy(hdr);
        sam_close(fp);
        return 1;
    }

    bam1_t *b = bam_init1();
    int ret;

    // BAMファイルからレコードを読み込み、Pileupに追加
    while ((ret = sam_read1(fp, hdr, b)) >= 0) {
        if (bam_plp_push(plp, b) < 0) {
            fprintf(stderr, "Error pushing to pileup.\n");
            break;
        }

        const bam_pileup1_t *p;
        int tid, n_plp;
        hts_pos_t pos;

        // Pileupエントリを取得して表示
        while ((p = bam_plp64_next(plp, &tid, &pos, &n_plp)) != NULL) {
            printf("Position %ld on %s: %d reads\n", pos, hdr->target_name[tid], n_plp);
            for (int i = 0; i < n_plp; i++) {
                printf("  Read: %s at QPOS: %d\n", bam_get_qname(p[i].b), p[i].qpos);
            }
        }
    }

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

    return 0;
}

t2.c

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

// Structure to store the BAM file and its header
typedef struct {
    samFile *fp;      // Pointer to the BAM file
    bam_hdr_t *hdr;   // Pointer to the BAM header
} sam_handle;

// Callback function: Reads a BAM record from the file
static int read_bam_record(void *data, bam1_t *b) {
    sam_handle *handle = (sam_handle *)data;
    return sam_read1(handle->fp, handle->hdr, b);
}

// Function to initialize and run pileup processing
void process_pileup(sam_handle *handle) {
    bam_plp_t plp = bam_plp_init(read_bam_record, handle);
    if (plp == NULL) {
        fprintf(stderr, "Failed to initialize pileup.\n");
        return;
    }

    const bam_pileup1_t *p;
    int tid, n_plp;
    hts_pos_t pos;

    // Perform pileup processing and print results
    while ((p = bam_plp_auto(plp, &tid, (int *)&pos, &n_plp)) != NULL) {
        printf("Position %d on %s: %d reads\n", (int)pos, 
               handle->hdr->target_name[tid], n_plp);
        for (int i = 0; i < n_plp; i++) {
            printf("  Read: %s at QPOS: %d\n", bam_get_qname(p[i].b), p[i].qpos);
        }
    }

    bam_plp_destroy(plp);  // Destroy the pileup iterator
}

// Main function to open the BAM file and run pileup
int main(int argc, char *argv[]) {
    if (argc < 2) {
        fprintf(stderr, "Usage: %s <bam_file>\n", argv[0]);
        return 1;
    }

    // Open the BAM file in binary read mode
    samFile *fp = sam_open(argv[1], "rb");
    if (fp == NULL) {
        fprintf(stderr, "Failed to open BAM file: %s\n", argv[1]);
        return 1;
    }

    // Read the BAM header
    bam_hdr_t *hdr = sam_hdr_read(fp);
    if (hdr == NULL) {
        fprintf(stderr, "Failed to read BAM header.\n");
        sam_close(fp);
        return 1;
    }

    // Create a handle for the BAM file and header
    sam_handle handle = {fp, hdr};

    // Process the pileup
    process_pileup(&handle);

    // Clean up resources
    bam_hdr_destroy(hdr);
    sam_close(fp);

    return 0;
}