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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
// Copyright ©2011-2012 The bíogo Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package pals
import (
"errors"
"fmt"
"github.com/biogo/biogo/alphabet"
"github.com/biogo/biogo/seq"
"github.com/biogo/biogo/seq/linear"
"github.com/biogo/biogo/util"
)
var (
binSize int = 1 << 10
minPadding int = 50
)
// A Packer collects a set of sequence into a Packed sequence.
type Packer struct {
packed *Packed
lastPad int
length int
}
type Packed struct {
*linear.Seq
seqMap
}
// Convert coordinates in a packed sequence into a feat.Feature.
func (pa *Packed) feature(from, to int, comp bool) (*Feature, error) {
if comp {
from, to = pa.Len()-to, pa.Len()-from
}
if from >= to {
return nil, errors.New("pals: from > to")
}
// DPHit coordinates sometimes over/underflow.
// This is a lazy hack to work around it, should really figure
// out what is going on.
if from < 0 {
from = 0
}
if to > pa.Len() {
to = pa.Len()
}
// Take midpoint of segment -- lazy hack again, endpoints
// sometimes under / overflow
bin := (from + to) / (2 * binSize)
binCount := (pa.Len() + binSize - 1) / binSize
if bin < 0 || bin >= binCount {
return nil, fmt.Errorf("pals: bin %d out of range 0..%d", bin, binCount-1)
}
contigIndex := pa.seqMap.binMap[bin]
if contigIndex < 0 || contigIndex >= len(pa.seqMap.contigs) {
return nil, fmt.Errorf("pals: contig %s index %d out of range 0..%d", pa.ID, contigIndex, len(pa.seqMap.contigs))
}
length := to - from
if length < 0 {
return nil, errors.New("pals: length < 0")
}
contig := pa.seqMap.contigs[contigIndex]
contigFrom := from - contig.from
contigTo := contigFrom + length
if contigFrom < 0 {
contigFrom = 0
}
if contigTo > contig.Len() {
contigTo = contig.Len()
}
return &Feature{
ID: contig.ID,
From: contigFrom,
To: contigTo,
Loc: Contig(contig.ID),
}, nil
}
// Create a new Packer.
func NewPacker(id string) *Packer {
return &Packer{
packed: &Packed{
Seq: &linear.Seq{Annotation: seq.Annotation{ID: id}},
seqMap: seqMap{},
},
}
}
// Pack a sequence into the Packed sequence. Returns a string giving diagnostic information.
func (pa *Packer) Pack(seq *linear.Seq) (string, error) {
if pa.packed.Alpha == nil {
pa.packed.Alpha = seq.Alpha
} else if pa.packed.Alpha != seq.Alpha {
return "", errors.New("pals: alphabet mismatch")
}
c := contig{Seq: seq}
padding := binSize - seq.Len()%binSize
if padding < minPadding {
padding += binSize
}
pa.length += pa.lastPad
c.from = pa.length
pa.length += seq.Len()
pa.lastPad = padding
m := &pa.packed.seqMap
bins := make([]int, (padding+seq.Len())/binSize)
for i := 0; i < len(bins); i++ {
bins[i] = len(m.contigs)
}
m.binMap = append(m.binMap, bins...)
m.contigs = append(m.contigs, c)
return fmt.Sprintf("%20s\t%10d\t%7d-%-d", seq.ID[:util.Min(20, len(seq.ID))], seq.Len(), len(m.binMap)-len(bins), len(m.binMap)-1), nil
}
// Finalise the sequence packing.
func (pa *Packer) FinalisePack() *Packed {
lastPad := 0
seq := make(alphabet.Letters, 0, pa.length)
for _, c := range pa.packed.seqMap.contigs {
padding := binSize - c.Len()%binSize
if padding < minPadding {
padding += binSize
}
seq = append(seq, alphabet.Letter('N').Repeat(lastPad)...)
seq = append(seq, c.Seq.Seq...)
lastPad = padding
}
pa.packed.Seq.Seq = seq
return pa.packed
}
// A contig holds a sequence within a SeqMap.
type contig struct {
*linear.Seq
from int
}
// A seqMap is a collection of sequences mapped to a Packed sequence.
type seqMap struct {
contigs []contig
binMap []int
}
|
