C++ sam_close函数代码示例

/******************************************************************************
*
*   quit, while performing some cleanup
*
*   int FLAG: What to free/close/etc.
*             0x1 things created by create_fastq_names()
*             0x2 things pthreads are closed and bam headers destroyed
*             In addition, the master node will free chromosomes.genome, close
*             the BAM file, and free everything in the chromosomes struct.
*
*   int rv: return value
*
*******************************************************************************/
void quit(int FLAG, int rv) {
    int taskid, i;

    free(config.bowtie2_options);

    MPI_Comm_rank(MPI_COMM_WORLD, &taskid);

    if(FLAG & 1) { //FASTQ filenames set
#ifndef DEBUG
        if(taskid == MASTER) {
            if(config.FASTQ1CT != NULL) remove(config.FASTQ1CT);
            if(config.paired && (config.FASTQ2GA != NULL)) remove(config.FASTQ2GA);
            if(!config.directional) {
                if(config.FASTQ1GA != NULL) remove(config.FASTQ1GA);
                if(config.paired && (config.FASTQ2CT != NULL)) remove(config.FASTQ2CT);
            }
        }
#endif
        if(config.FASTQ1CT != NULL) free(config.FASTQ1CT);
        if(config.FASTQ1GA != NULL) free(config.FASTQ1GA);
        if(config.unmapped1 != NULL) free(config.unmapped1);
        if(config.paired) {
            if(config.FASTQ2CT != NULL) free(config.FASTQ2CT);
            if(config.FASTQ2GA != NULL) free(config.FASTQ2GA);
            if(config.unmapped2 != NULL) free(config.unmapped2);
        }
        free(config.basename);
        free(config.outname);
        if(config.fai) free(config.fai);
    }

    if(taskid == MASTER) {
        free(chromosomes.genome);
        for(i=0; i<chromosomes.nchromosomes; i++) {
            free((chromosomes.chromosome[i])->chrom);
            free(*(chromosomes.chromosome+i));
        }
        free(chromosomes.chromosome);
        if(FLAG && OUTPUT_BAM) sam_close(OUTPUT_BAM);
    }
    MPI_Finalize();
    if(taskid == MASTER && FLAG > 0) {
#ifdef DEBUG
        if(fp1) sam_close(fp1);
        if(fp2) sam_close(fp2);
        if(!config.directional) {
            if(fp3) sam_close(fp3);
            if(fp4) sam_close(fp4);
        }
#else
        if(config.unmapped) {
            pclose(unmapped1);
            if(config.paired) pclose(unmapped2);
        }
#endif
    }
    exit(rv);
}

static void cleanup_state(state_t* state)
{
    if (!state) return;
    free(state->rg_id);
    if (state->output_file) sam_close(state->output_file);
    bam_hdr_destroy(state->output_header);
    if (state->input_file) sam_close(state->input_file);
    bam_hdr_destroy(state->input_header);
    free(state);
}

static void copy_check_alignment(const char *infname, const char *informat,
    const char *outfname, const char *outmode, const char *outref)
{
    samFile *in = sam_open(infname, "r");
    samFile *out = sam_open(outfname, outmode);
    bam1_t *aln = bam_init1();
    bam_hdr_t *header = NULL;
    int res;

    if (!in) {
        fail("couldn't open %s", infname);
        goto err;
    }
    if (!out) {
        fail("couldn't open %s with mode %s", outfname, outmode);
        goto err;
    }
    if (!aln) {
        fail("bam_init1() failed");
        goto err;
    }

    if (outref) {
        if (hts_set_opt(out, CRAM_OPT_REFERENCE, outref) < 0) {
            fail("setting reference %s for %s", outref, outfname);
            goto err;
        }
    }

    header = sam_hdr_read(in);
    if (!header) {
        fail("reading header from %s", infname);
        goto err;
    }
    if (sam_hdr_write(out, header) < 0) fail("writing headers to %s", outfname);

    while ((res = sam_read1(in, header, aln)) >= 0) {
        int mod4 = ((intptr_t) bam_get_cigar(aln)) % 4;
        if (mod4 != 0)
            fail("%s CIGAR not 4-byte aligned; offset is 4k+%d for \"%s\"",
                 informat, mod4, bam_get_qname(aln));

        if (sam_write1(out, header, aln) < 0) fail("writing to %s", outfname);
    }
    if (res < -1) {
        fail("failed to read alignment from %s", infname);
    }

 err:
    bam_destroy1(aln);
    bam_hdr_destroy(header);
    if (in) sam_close(in);
    if (out) sam_close(out);
}

samfile_t *samopen(const char *fn, const char *mode, const void *aux)
{
    // hts_open() is really sam_open(), except for #define games
    samFile *hts_fp = hts_open(fn, mode);
    if (hts_fp == NULL)  return NULL;

    samfile_t *fp = malloc(sizeof (samfile_t));
    if (!fp) {
        sam_close(hts_fp);
        return NULL;
    }
    fp->file = hts_fp;
    fp->x.bam = hts_fp->fp.bgzf;
    if (strchr(mode, 'r')) {
        if (aux) {
            if (hts_set_fai_filename(fp->file, aux) != 0) {
                sam_close(hts_fp);
                free(fp);
                return NULL;
            }
        }
        fp->header = sam_hdr_read(fp->file);  // samclose() will free this
        if (fp->header == NULL) {
            sam_close(hts_fp);
            free(fp);
            return NULL;
        }
        fp->is_write = 0;
        if (fp->header->n_targets == 0 && bam_verbose >= 1)
            fprintf(samtools_stderr, "[samopen] no @SQ lines in the header.\n");
    }
    else {
        enum htsExactFormat fmt = hts_get_format(fp->file)->format;
        fp->header = (bam_hdr_t *)aux;  // For writing, we won't free it
        fp->is_write = 1;
        if (!(fmt == text_format || fmt == sam) || strchr(mode, 'h')) {
            if (sam_hdr_write(fp->file, fp->header) < 0) {
                if (bam_verbose >= 1)
                    fprintf(samtools_stderr, "[samopen] Couldn't write header\n");
                sam_close(hts_fp);
                free(fp);
                return NULL;
            }
        }
    }

    return fp;
}

int main(int argc, char **argv) {
    dlib::BamHandle in = dlib::BamHandle("bed_test.bam");
    dlib::ParsedBed bed = dlib::ParsedBed("bed_test.bed", in.header);
    bam1_t *b = bam_init1();
    size_t diffs = 0;
    void *lh3bed = bed_read("bed_test.bed");
    samFile *so = sam_open("disagreed.bam", "wb9");
    sam_hdr_write(so, in.header);
    size_t disagrees = 0, agrees = 0;
    int dbr = 0, lh3r = 0;
    while(in.read(b) != -1) {
        if(b->core.flag & (BAM_FUNMAP)) continue;
        if((dbr = bed.bam1_test(b)) != (lh3r = bed_overlap(lh3bed, in.header->target_name[b->core.tid], b->core.pos, bam_endpos(b)))) {
            LOG_EXIT("dbr: %i. lh3r: %i. Contig: %s. Position: %i. endpos; %i\n", dbr, lh3r, in.header->target_name[b->core.tid], b->core.pos, bam_endpos(b));
            if(++disagrees % 100 == 0) LOG_DEBUG("disagrees: %lu.\n", disagrees);
            sam_write1(so, in.header, b);
        } else {
            if(++agrees % 500000 == 0) LOG_DEBUG("agrees: %lu.\n", agrees);
        }
    }
    sam_close(so);
    bam_destroy1(b);
    bed_destroy(lh3bed);
    return EXIT_SUCCESS;
}

int main_cat(int argc, char *argv[])
{
    bam_header_t *h = 0;
    char *outfn = 0;
    int c, ret;
    while ((c = getopt(argc, argv, "h:o:")) >= 0) {
        switch (c) {
            case 'h': {
                tamFile fph = sam_open(optarg);
                if (fph == 0) {
                    fprintf(stderr, "[%s] ERROR: fail to read the header from '%s'.\n", __func__, argv[1]);
                    return 1;
                }
                h = sam_header_read(fph);
                sam_close(fph);
                break;
            }
            case 'o': outfn = strdup(optarg); break;
        }
    }
    if (argc - optind < 2) {
        fprintf(stderr, "Usage: samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [...]\n");
        return 1;
    }
    ret = bam_cat(argc - optind, argv + optind, h, outfn? outfn : "-");
    free(outfn);
    return ret;
}

int main_reheader(int argc, char *argv[])
{
	bam_header_t *h;
	BGZF *in;
	if (argc != 3) {
		fprintf(stderr, "Usage: samtools reheader <in.header.sam> <in.bam>\n");
		return 1;
	}
	{ // read the header
		tamFile fph = sam_open(argv[1]);
		if (fph == 0) {
			fprintf(stderr, "[%s] fail to read the header from %s.\n", __func__, argv[1]);
			return 1;
		}
		h = sam_header_read(fph);
		sam_close(fph);
	}
	in = strcmp(argv[2], "-")? bam_open(argv[2], "r") : bam_dopen(fileno(stdin), "r");
	if (in == 0) {
		fprintf(stderr, "[%s] fail to open file %s.\n", __func__, argv[2]);
		return 1;
	}
	bam_reheader(in, h, fileno(stdout));
	bgzf_close(in);
	return 0;
}

/*static*/ bam_header_t * SAM::update_header_from_list(bam_header_t *header, names_list_t & list) {
	Temporary_File samfile;
	samfile.close_file();
	samfile_t * sf = samopen(samfile.get_filename().c_str(),"wh",header);
	samclose(sf);

	Temporary_File tempfile;
	ofstream &output = tempfile.get_stream();

	ifstream input(samfile.get_filename().c_str());
	string temp;
	while (not input.eof()) {
		getline(input,temp);
		if ((temp.size() >= 3) and (temp[0] != '@' or temp[1] != 'S' or temp[2] != 'Q'))
			output << temp << '\n';
	}

	for (names_list_t::iterator iter = list.begin(); iter != list.end(); iter++)
		output << "@SQ\tSN:" << iter->first << "\tLN:" << iter->second << '\n';
	tempfile.close_file();

	tamFile fp = sam_open(tempfile.get_filename().c_str());

	bam_header_t * newheader = sam_header_read(fp);
	sam_close(fp);

	return newheader;
}

void metaBigClose(struct metaBig** pMb)
/* close the file and free up everything. */
{
    struct metaBig* mb = *pMb;
    hashFree(&mb->chromSizeHash);
    if (mb->rgList)
        hashFree(&mb->rgList);
    if (mb->sections)
        bedFreeList(&mb->sections);
    if (mb->originalFileName)
        freeMem(mb->originalFileName);
    if (mb->fileName)
        freeMem(mb->fileName);
    if (mb->baseFileName)
        freeMem(mb->baseFileName);
    if (mb->remoteSiteAndDir)
        freeMem(mb->remoteSiteAndDir);
#ifdef USE_HTSLIB
    if (mb->idx)
        hts_idx_destroy(mb->idx);
#endif
    if (mb->type == isaBigBed)
        bigBedFileClose(&mb->big.bbi);
#ifdef USE_HTSLIB
    else if (mb->type == isaBam)
        sam_close(mb->big.bam);
#endif
    else
        bigWigFileClose(&mb->big.bbi);
#ifdef USE_HTSLIB
    if (mb->header)
        bam_hdr_destroy(mb->header);
#endif
    freez(pMb);
}

int bam_idxstats(int argc, char *argv[])
{
    hts_idx_t* idx;
    bam_hdr_t* header;
    samFile* fp;

    if (argc < 2) {
        fprintf(pysamerr, "Usage: samtools idxstats <in.bam>\n");
        return 1;
    }
    fp = sam_open(argv[1], "r");
    if (fp == NULL) { fprintf(pysamerr, "[%s] fail to open BAM.\n", __func__); return 1; }
    header = sam_hdr_read(fp);
    idx = sam_index_load(fp, argv[1]);
    if (idx == NULL) { fprintf(pysamerr, "[%s] fail to load the index.\n", __func__); return 1; }

    int i;
    for (i = 0; i < header->n_targets; ++i) {
        // Print out contig name and length
        printf("%s\t%d", header->target_name[i], header->target_len[i]);
        // Now fetch info about it from the meta bin
        uint64_t u, v;
        hts_idx_get_stat(idx, i, &u, &v);
        printf("\t%" PRIu64 "\t%" PRIu64 "\n", u, v);
    }
    // Dump information about unmapped reads
    printf("*\t0\t0\t%" PRIu64 "\n", hts_idx_get_n_no_coor(idx));
    bam_hdr_destroy(header);
    hts_idx_destroy(idx);
    sam_close(fp);
    return 0;
}

// remember to clean up with bam_destroy1(b);
bam1_t* alignment_to_bam(const string& sam_header,
                         const Alignment& alignment,
                         const string& refseq,
                         const int32_t refpos,
                         const string& cigar,
                         const string& mateseq,
                         const int32_t matepos,
                         const int32_t tlen) {

    assert(!sam_header.empty());
    string sam_file = "data:" + sam_header + alignment_to_sam(alignment, refseq, refpos, cigar, mateseq, matepos, tlen);
    const char* sam = sam_file.c_str();
    samFile *in = sam_open(sam, "r");
    bam_hdr_t *header = sam_hdr_read(in);
    bam1_t *aln = bam_init1();
    if (sam_read1(in, header, aln) >= 0) {
        bam_hdr_destroy(header);
        sam_close(in); // clean up
        return aln;
    } else {
        cerr << "[vg::alignment] Failure to parse SAM record" << endl
             << sam << endl;
        exit(1);
    }
}

static int aux_fields1(void)
{
    static const char sam[] = "data:"
"@SQ\tSN:one\tLN:1000\n"
"@SQ\tSN:two\tLN:500\n"
"r1\t0\tone\t500\t20\t8M\t*\t0\t0\tATGCATGC\tqqqqqqqq\tXA:A:k\tXi:i:37\tXf:f:" xstr(PI) "\tXd:d:" xstr(E) "\tXZ:Z:" HELLO "\tXH:H:" BEEF "\tXB:B:c,-2,0,+2\tZZ:i:1000000\n";

    // Canonical form of the alignment record above, as output by sam_format1()
    static const char r1[] = "r1\t0\tone\t500\t20\t8M\t*\t0\t0\tATGCATGC\tqqqqqqqq\tXA:A:k\tXi:i:37\tXf:f:3.14159\tXd:d:2.71828\tXZ:Z:" HELLO "\tXH:H:" BEEF "\tXB:B:c,-2,0,2\tZZ:i:1000000";

    samFile *in = sam_open(sam, "r");
    bam_hdr_t *header = sam_hdr_read(in);
    bam1_t *aln = bam_init1();
    uint8_t *p;
    uint32_t n;
    kstring_t ks = { 0, 0, NULL };

    if (sam_read1(in, header, aln) >= 0) {
        if ((p = check_bam_aux_get(aln, "XA", 'A')) && bam_aux2A(p) != 'k')
            fail("XA field is '%c', expected 'k'", bam_aux2A(p));

        if ((p = check_bam_aux_get(aln, "Xi", 'C')) && bam_aux2i(p) != 37)
            fail("Xi field is %d, expected 37", bam_aux2i(p));

        if ((p = check_bam_aux_get(aln, "Xf", 'f')) && fabs(bam_aux2f(p) - PI) > 1E-6)
            fail("Xf field is %.12f, expected pi", bam_aux2f(p));

        if ((p = check_bam_aux_get(aln, "Xd", 'd')) && fabs(bam_aux2f(p) - E) > 1E-6)
            fail("Xf field is %.12f, expected e", bam_aux2f(p));

        if ((p = check_bam_aux_get(aln, "XZ", 'Z')) && strcmp(bam_aux2Z(p), HELLO) != 0)
            fail("XZ field is \"%s\", expected \"%s\"", bam_aux2Z(p), HELLO);

        if ((p = check_bam_aux_get(aln, "XH", 'H')) && strcmp(bam_aux2Z(p), BEEF) != 0)
            fail("XH field is \"%s\", expected \"%s\"", bam_aux2Z(p), BEEF);

        // TODO Invent and use bam_aux2B()
        if ((p = check_bam_aux_get(aln, "XB", 'B')) && ! (memcmp(p, "Bc", 2) == 0 && (memcpy(&n, p+2, 4), n) == 3 && memcmp(p+6, "\xfe\x00\x02", 3) == 0))
            fail("XB field is %c,..., expected c,-2,0,+2", p[1]);

        if ((p = check_bam_aux_get(aln, "ZZ", 'I')) && bam_aux2i(p) != 1000000)
            fail("ZZ field is %d, expected 1000000", bam_aux2i(p));

        if (sam_format1(header, aln, &ks) < 0)
            fail("can't format record");

        if (strcmp(ks.s, r1) != 0)
            fail("record formatted incorrectly: \"%s\"", ks.s);

        free(ks.s);
    }
    else fail("can't read record");

    bam_destroy1(aln);
    bam_hdr_destroy(header);
    sam_close(in);

    return 1;
}

void samclose(samfile_t *fp)
{
    if (fp == 0) return;
    if (fp->header) bam_header_destroy(fp->header);
    if (fp->type & 1) bam_close(fp->x.bam);
    else if (fp->type == 2) sam_close(fp->x.tamr);
    free(fp);
}

void samclose(samfile_t *fp)
{
    if (fp) {
        if (!fp->is_write && fp->header) bam_hdr_destroy(fp->header);
        sam_close(fp->file);
        free(fp);
    }
}

int bam_mating(int argc, char *argv[])
{
    samFile *in, *out;
    int c, remove_reads = 0, proper_pair_check = 1, add_ct = 0;
    sam_global_args ga = SAM_GLOBAL_ARGS_INIT;
    char wmode[3] = {'w', 'b', 0};
    static const struct option lopts[] = {
        SAM_OPT_GLOBAL_OPTIONS('-', 0, 'O', 0, 0),
        { NULL, 0, NULL, 0 }
    };

    // parse args
    if (argc == 1) { usage(stdout); return 0; }
    while ((c = getopt_long(argc, argv, "rpcO:", lopts, NULL)) >= 0) {
        switch (c) {
            case 'r': remove_reads = 1; break;
            case 'p': proper_pair_check = 0; break;
            case 'c': add_ct = 1; break;
            default:  if (parse_sam_global_opt(c, optarg, lopts, &ga) == 0) break;
                      /* else fall-through */
            case '?': usage(stderr); return 1;
        }
    }
    if (optind+1 >= argc) { usage(stderr); return 1; }

    // init
    if ((in = sam_open_format(argv[optind], "rb", &ga.in)) == NULL) {
        fprintf(stderr, "[bam_mating] cannot open input file\n");
        return 1;
    }
    sam_open_mode(wmode+1, argv[optind+1], NULL);
    if ((out = sam_open_format(argv[optind+1], wmode, &ga.out)) == NULL) {
        fprintf(stderr, "[bam_mating] cannot open output file\n");
        return 1;
    }

    // run
    bam_mating_core(in, out, remove_reads, proper_pair_check, add_ct);

    // cleanup
    sam_close(in); sam_close(out);
    sam_global_args_free(&ga);

    return 0;
}