forked from samtools/samtools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sam.c
83 lines (76 loc) · 2.26 KB
/
sam.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include <string.h>
#include <unistd.h>
#include "htslib/faidx.h"
#include "sam.h"
int samthreads(samfile_t *fp, int n_threads, int n_sub_blks)
{
if (!fp->file->is_bin || !fp->file->is_write) return -1;
bgzf_mt(fp->x.bam, n_threads, n_sub_blks);
return 0;
}
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));
fp->file = hts_fp;
fp->x.bam = hts_fp->fp.bgzf;
if (strchr(mode, 'r')) {
if (aux) hts_set_fai_filename(fp->file, aux);
fp->header = sam_hdr_read(fp->file); // samclose() will free this
if (fp->header->n_targets == 0 && bam_verbose >= 1)
fprintf(stderr, "[samopen] no @SQ lines in the header.\n");
}
else {
fp->header = (bam_hdr_t *)aux; // For writing, we won't free it
if (fp->file->is_bin || fp->file->is_cram || strchr(mode, 'h')) sam_hdr_write(fp->file, fp->header);
}
return fp;
}
void samclose(samfile_t *fp)
{
if (fp) {
if (!fp->file->is_write && fp->header) bam_hdr_destroy(fp->header);
sam_close(fp->file);
free(fp);
}
}
int sampileup(samfile_t *fp, int mask, bam_pileup_f func, void *func_data)
{
bam_plbuf_t *buf;
int ret;
bam1_t *b;
b = bam_init1();
buf = bam_plbuf_init(func, func_data);
if (mask < 0) mask = BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP;
else mask |= BAM_FUNMAP;
while ((ret = samread(fp, b)) >= 0) {
// bam_plp_push() itself now filters out unmapped reads only
if (b->core.flag & mask) b->core.flag |= BAM_FUNMAP;
bam_plbuf_push(b, buf);
}
bam_plbuf_push(0, buf);
bam_plbuf_destroy(buf);
bam_destroy1(b);
return 0;
}
char *samfaipath(const char *fn_ref)
{
char *fn_list = 0;
if (fn_ref == 0) return 0;
fn_list = calloc(strlen(fn_ref) + 5, 1);
strcat(strcpy(fn_list, fn_ref), ".fai");
if (access(fn_list, R_OK) == -1) { // fn_list is unreadable
if (access(fn_ref, R_OK) == -1) {
fprintf(stderr, "[samfaipath] fail to read file %s.\n", fn_ref);
} else {
if (bam_verbose >= 3) fprintf(stderr, "[samfaipath] build FASTA index...\n");
if (fai_build(fn_ref) == -1) {
fprintf(stderr, "[samfaipath] fail to build FASTA index.\n");
free(fn_list); fn_list = 0;
}
}
}
return fn_list;
}