Files
obitools4/pkg/obiapat/pcr.go
2022-01-14 17:17:54 +01:00

428 lines
10 KiB
Go

package obiapat
import (
"git.metabarcoding.org/lecasofts/go/obitools/pkg/goutils"
"git.metabarcoding.org/lecasofts/go/obitools/pkg/obiseq"
)
type _Options struct {
minLength int
maxLength int
circular bool
forwardError int
reverseError int
bufferSize int
batchSize int
parallelWorkers int
}
// Options stores a set of option usable by the
// PCR simulation algotithm.
type Options struct {
pointer *_Options
}
// WithOption is the standard type for function
// declaring options.
type WithOption func(Options)
// MinLength method returns minimum length of
// the searched amplicon (length of the primers
// excluded)
func (options Options) MinLength() int {
return options.pointer.minLength
}
// MaxLength method returns maximum length of
// the searched amplicon (length of the primers
// excluded)
func (options Options) MaxLength() int {
return options.pointer.maxLength
}
// ForwardError method returns the number of
// error allowed when matching the forward
// primer.
func (options Options) ForwardError() int {
return options.pointer.forwardError
}
// ReverseError method returns the number of
// error allowed when matching the reverse
// primer.
func (options Options) ReverseError() int {
return options.pointer.reverseError
}
// Circular method returns the topology option.
// true for circular, false for linear
func (options Options) Circular() bool {
return options.pointer.circular
}
// BufferSize returns the size of the channel
// buffer specified by the options
func (options Options) BufferSize() int {
return options.pointer.bufferSize
}
// BatchSize returns the size of the
// sequence batch used by the PCR algorithm
func (options Options) BatchSize() int {
return options.pointer.batchSize
}
// ParallelWorkers returns how many search
// jobs will be run in parallel.
func (options Options) ParallelWorkers() int {
return options.pointer.parallelWorkers
}
// MakeOptions buils a new default option set for
// the PCR simulation algoithm.
func MakeOptions(setters []WithOption) Options {
o := _Options{
minLength: 0,
maxLength: 0,
forwardError: 0,
reverseError: 0,
circular: false,
parallelWorkers: 4,
batchSize: 100,
bufferSize: 100,
}
opt := Options{&o}
for _, set := range setters {
set(opt)
}
return opt
}
// OptionMinLength sets the minimum length of
// the searched amplicon (length of the primers
// excluded)
func OptionMinLength(minLength int) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.minLength = minLength
})
return f
}
// OptionMaxLength sets the maximum length of
// the searched amplicon (length of the primers
// excluded)
func OptionMaxLength(maxLength int) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.maxLength = maxLength
})
return f
}
// OptionForwardError sets the number of
// error allowed when matching the forward
// primer.
func OptionForwardError(max int) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.forwardError = max
})
return f
}
// OptionReverseError sets the number of
// error allowed when matching the reverse
// primer.
func OptionReverseError(max int) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.reverseError = max
})
return f
}
// OptionCircular sets the topology option.
// true for circular, false for linear
func OptionCircular(circular bool) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.circular = circular
})
return f
}
// OptionBufferSize sets the requested channel
// buffer size.
func OptionBufferSize(size int) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.bufferSize = size
})
return f
}
// OptionParallelWorkers sets how many search
// jobs will be run in parallel.
func OptionParallelWorkers(nworkers int) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.parallelWorkers = nworkers
})
return f
}
// OptionBatchSize sets the requested sequence
// batch size.
func OptionBatchSize(size int) WithOption {
f := WithOption(func(opt Options) {
opt.pointer.batchSize = size
})
return f
}
func _Pcr(seq ApatSequence, sequence obiseq.BioSequence,
forward, cfwd, reverse, crev ApatPattern,
opt Options) obiseq.BioSequenceSlice {
results := make(obiseq.BioSequenceSlice, 0, 10)
forwardMatches := forward.FindAllIndex(seq)
if forwardMatches != nil {
begin := forwardMatches[0][0]
length := seq.Length() - begin
if opt.pointer.maxLength > 0 {
length = forwardMatches[len(forwardMatches)-1][2] - begin + opt.MaxLength() + reverse.Length()
}
if opt.Circular() {
begin = 0
length = seq.Length() + _MaxPatLen
}
reverseMatches := crev.FindAllIndex(seq, begin, length)
if reverseMatches != nil {
for _, fm := range forwardMatches {
posi := fm[0]
if posi < seq.Length() {
erri := fm[2]
for _, rm := range reverseMatches {
posj := rm[0]
if posj < seq.Length() {
posj := rm[1]
errj := rm[2]
length = 0
if posj > posi {
length = rm[0] - fm[1]
} else {
if opt.Circular() {
length = rm[0] + seq.Length() - posi - forward.Length()
}
}
if length > 0 && // For when primers touch or overlap
(opt.MinLength() == 0 || length >= opt.MinLength()) &&
(opt.MaxLength() == 0 || length <= opt.MaxLength()) {
amplicon, _ := sequence.Subsequence(fm[1], rm[0], opt.pointer.circular)
annot := amplicon.Annotations()
goutils.CopyMap(annot, sequence.Annotations())
annot["forward_primer"] = forward.String()
match, _ := sequence.Subsequence(fm[0], fm[1], opt.pointer.circular)
annot["forward_match"] = match.String()
match.Destroy()
annot["forward_error"] = erri
annot["reverse_primer"] = reverse.String()
match, _ = sequence.Subsequence(rm[0], rm[1], opt.pointer.circular)
match = match.ReverseComplement(true)
annot["reverse_match"] = match.String()
match.Destroy()
annot["reverse_error"] = errj
results = append(results, amplicon)
}
}
}
}
}
}
}
forwardMatches = reverse.FindAllIndex(seq)
if forwardMatches != nil {
begin := forwardMatches[0][0]
length := seq.Length() - begin
if opt.pointer.maxLength > 0 {
length = forwardMatches[len(forwardMatches)-1][2] - begin + opt.MaxLength() + reverse.Length()
}
if opt.Circular() {
begin = 0
length = seq.Length() + _MaxPatLen
}
reverseMatches := cfwd.FindAllIndex(seq, begin, length)
if reverseMatches != nil {
for _, fm := range forwardMatches {
posi := fm[0]
if posi < seq.Length() {
erri := fm[2]
for _, rm := range reverseMatches {
posj := rm[0]
if posj < seq.Length() {
posj := rm[1]
errj := rm[2]
length = 0
if posj > posi {
length = rm[0] - fm[1]
} else {
if opt.Circular() {
length = rm[0] + seq.Length() - posi - forward.Length()
}
}
if length > 0 && // For when primers touch or overlap
(opt.MinLength() == 0 || length >= opt.MinLength()) &&
(opt.MaxLength() == 0 || length <= opt.MaxLength()) {
amplicon, _ := sequence.Subsequence(fm[1], rm[0], opt.pointer.circular)
amplicon = amplicon.ReverseComplement(true)
annot := amplicon.Annotations()
goutils.CopyMap(annot, sequence.Annotations())
annot["forward_primer"] = forward.String()
match, _ := sequence.Subsequence(rm[0], rm[1], opt.pointer.circular)
match.ReverseComplement(true)
annot["forward_match"] = match.String()
match.Destroy()
annot["forward_error"] = errj
annot["reverse_primer"] = reverse.String()
match, _ = sequence.Subsequence(fm[0], fm[1], opt.pointer.circular)
annot["reverse_match"] = match.String()
match.Destroy()
annot["reverse_error"] = erri
results = append(results, amplicon)
}
}
}
}
}
}
}
return results
}
// PCR runs the PCR simulation algorithm on a single
// obiseq.BioSequence instance. PCR parameters are
// specified using the corresponding Option functions
// defined for the PCR algorithm.
func PCR(sequence obiseq.BioSequence,
forward, reverse string, options ...WithOption) obiseq.BioSequenceSlice {
opt := MakeOptions(options)
seq, _ := MakeApatSequence(sequence, opt.Circular())
fwd, _ := MakeApatPattern(forward, opt.ForwardError())
rev, _ := MakeApatPattern(reverse, opt.ReverseError())
cfwd, _ := fwd.ReverseComplement()
crev, _ := rev.ReverseComplement()
results := _Pcr(seq, sequence,
fwd, cfwd, rev, crev,
opt)
seq.Free()
fwd.Free()
rev.Free()
cfwd.Free()
crev.Free()
return results
}
// PCRSlice runs the PCR simulation algorithm on a set of
// obiseq.BioSequence instances grouped in a obiseq.BioSequenceSlice.
// PCR parameters are
// specified using the corresponding Option functions
// defined for the PCR algorithm.
func PCRSlice(sequences obiseq.BioSequenceSlice,
forward, reverse string, options ...WithOption) obiseq.BioSequenceSlice {
results := make(obiseq.BioSequenceSlice, 0, len(sequences))
opt := MakeOptions(options)
fwd, _ := MakeApatPattern(forward, opt.ForwardError())
rev, _ := MakeApatPattern(reverse, opt.ReverseError())
cfwd, _ := fwd.ReverseComplement()
crev, _ := rev.ReverseComplement()
if len(sequences) > 0 {
seq, _ := MakeApatSequence(sequences[0], opt.Circular())
amplicons := _Pcr(seq, sequences[0],
fwd, cfwd, rev, crev,
opt)
if len(amplicons) > 0 {
results = append(results, amplicons...)
}
for _, sequence := range sequences[1:] {
seq, _ := MakeApatSequence(sequence, opt.Circular(), seq)
amplicons = _Pcr(seq, sequence,
fwd, cfwd, rev, crev,
opt)
if len(amplicons) > 0 {
results = append(results, amplicons...)
}
}
seq.Free()
}
fwd.Free()
rev.Free()
cfwd.Free()
crev.Free()
return results
}
// PCRSliceWorker is a worker function builder which produce
// job function usable by the obiseq.MakeISliceWorker function.
func PCRSliceWorker(forward, reverse string,
options ...WithOption) obiseq.SeqSliceWorker {
worker := func(sequences obiseq.BioSequenceSlice) obiseq.BioSequenceSlice {
return PCRSlice(sequences, forward, reverse, options...)
}
return worker
}