forked from brentp/bigly
/
bamat.go
89 lines (81 loc) · 1.76 KB
/
bamat.go
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
84
85
86
87
88
89
package bamat
import (
"bufio"
"os"
"github.com/biogo/hts/bam"
"github.com/biogo/hts/sam"
)
// BamAt holds bam index and the Bam Reader.
// Because BamAt holds the underlying os.File open, it is not
// safe to query from multiple go routines.
type BamAt struct {
*bam.Reader
idx *bam.Index
fh *os.File
Refs map[string]*sam.Reference
}
// New returns a BamAt from the given path to the indexed bam
func New(path string) (*BamAt, error) {
bamat := &BamAt{}
if !(path == "" || path == "-" || path == "stdin") {
f, err := os.Open(path + ".bai")
if err != nil && len(path) > 4 {
f, err = os.Open(path[:len(path)-4] + ".bai")
}
if err != nil {
return nil, err
}
defer f.Close()
idx, err := bam.ReadIndex(bufio.NewReader(f))
if err != nil {
return nil, err
}
bamat.idx = idx
bamat.fh, err = os.Open(path)
if err != nil {
return nil, err
}
} else {
bamat.fh = os.Stdin
}
br, err := bam.NewReader(bamat.fh, 2)
if err != nil {
bamat.fh.Close()
return nil, err
}
bamat.Reader = br
hdr := br.Header()
bamat.Refs = make(map[string]*sam.Reference, 40)
for _, r := range hdr.Refs() {
bamat.Refs[r.Name()] = r
}
return bamat, nil
}
// Query the BamAt with 0-base half-open interval.
func (b *BamAt) Query(chrom string, start int, end int) (*bam.Iterator, error) {
if chrom == "" { // stdin
return bam.NewIterator(b.Reader, nil)
}
ref := b.Refs[chrom]
if end <= 0 {
end = ref.Len() - 1
}
chunks, err := b.idx.Chunks(ref, start, end)
if err != nil {
return nil, err
}
return bam.NewIterator(b.Reader, chunks)
}
// Close closes the underlying file and the Bam.Reader
func (b *BamAt) Close() error {
if b == nil {
return nil
}
if b.Reader != nil {
b.Reader.Close()
}
if b.fh != nil {
return b.fh.Close()
}
return nil
}