Browse Source

basesWith: MD tag

remotes/georgeg/no_streams
lomereiter 10 years ago
parent
commit
ec2c9aeb0c
  1. 175
      bio/bam/baseinfo.d
  2. 115
      bio/bam/md/operation.d
  3. 139
      bio/bam/md/parse.d
  4. 2
      bio/bam/md/reconstruct.d
  5. 30
      bio/core/base.d
  6. 223
      bio/core/sequence.d
  7. 24
      test/unittests.d

175
bio/bam/baseinfo.d

@ -20,9 +20,12 @@
module bio.bam.baseinfo;
import bio.core.base;
import bio.core.sequence;
import bio.bam.read;
import bio.bam.tagvalue;
import bio.bam.fz.flowcall;
import bio.bam.md.core;
import std.range;
import std.conv;
@ -67,6 +70,7 @@ struct PerBaseInfo(R, Tags...) {
// - void populate(Result)(ref Result r);
// Populates fields of the result declared in resultProperties.
// Should run in O(1), just copying a few variables.
// Current base of the result is updated before the call.
//
// - void update(const ref R read);
// Encapsulates logic of moving to the next base and updating
@ -86,6 +90,7 @@ struct PerBaseInfo(R, Tags...) {
}
static struct Result {
/// Actual read base, with strand taken into account.
Base base;
alias base this;
@ -124,21 +129,21 @@ struct PerBaseInfo(R, Tags...) {
this(Args...)(R read, Args args) {
_read = read;
_seq = read.sequence;
_rev = read.is_reverse_strand;
_seq = reversableRange!complementBase(read.sequence, _rev);
foreach (t; Extensions) {
setup!t(read, args);
}
}
bool empty() @property const {
bool empty() @property {
return _seq.empty;
}
Result front() @property {
Result r = void;
r.base = _rev ? _seq.back.complement : _seq.front;
r.base = _seq.front;
foreach (t; Extensions)
populate!t(r);
return r;
@ -151,7 +156,7 @@ struct PerBaseInfo(R, Tags...) {
PerBaseInfo!(R, Tags) save() @property {
PerBaseInfo!(R, Tags) r = void;
r._read = _read.dup;
r._seq = r._read.sequence;
r._seq = _seq.save;
r._rev = _rev;
foreach (t; Extensions)
copy!t(r);
@ -164,15 +169,18 @@ struct PerBaseInfo(R, Tags...) {
update!t();
}
if (_rev)
_seq.popBack();
else
_seq.popFront();
_seq.popFront();
}
/// Returns true if the read is reverse strand,
/// and false otherwise.
bool reverse_strand() @property const {
return _rev;
}
private {
R _read = void;
typeof(_read.sequence) _seq = void;
ReversableRange!(complementBase, typeof(_read.sequence)) _seq = void;
bool _rev = void;
}
}
@ -191,6 +199,108 @@ template basesWith(Tags...) {
}
}
/// Provides additional properties
/// * reference_base
/// * md_operation
template MDbaseInfo(R) {
mixin template resultProperties() {
/// If current CIGAR operation is reference consuming,
/// returns reference base at this position, otherwise
/// returns '-'.
///
/// If read is on '-' strand, the result will be
/// complementary base.
char reference_base() @property {
return _ref_base;
}
/// Current MD operation
MdOperation md_operation() @property {
return _md_op;
}
private {
char _ref_base = void;
MdOperation _md_op = void;
}
}
mixin template rangeMethods() {
private {
ReversableRange!(reverseMdOp, MdOperationRange) _md_ops = void;
uint _match; // remaining length of current match operation
}
void setup(Args...)(const ref R read, Args args)
{
auto md = read["MD"];
auto md_str = *(cast(string*)&md);
_md_ops = reversableRange!reverseMdOp(mdOperations(md_str),
read.is_reverse_strand);
while (!_md_ops.empty && _md_ops.front.is_deletion)
_md_ops.popFront();
if (!_md_ops.empty && _md_ops.front.is_match)
_match = _md_ops.front.match;
}
void populate(Result)(ref Result result)
{
if (!current_cigar_operation.is_reference_consuming)
{
result._ref_base = '-';
return;
}
if (_md_ops.empty) {
return;
}
auto op = _md_ops.front;
if (op.is_mismatch)
result._ref_base = op.mismatch.asCharacter;
else if (op.is_match) {
result._ref_base = result.base.asCharacter;
}
else assert(0);
result._md_op = op;
}
void update(const ref R read)
{
if (!current_cigar_operation.is_reference_consuming)
return;
if (_md_ops.front.is_mismatch)
{
_md_ops.popFront();
}
else if (_md_ops.front.is_match)
{
--_match;
if (_match == 0)
_md_ops.popFront();
}
else assert(0);
while (!_md_ops.empty && _md_ops.front.is_deletion)
_md_ops.popFront();
if (_match == 0 && !_md_ops.empty && _md_ops.front.is_match)
_match = _md_ops.front.match;
}
void copy(Range)(ref Range source, ref Range target)
{
target.MD._md_ops = source.MD._md_ops.save;
}
}
}
/// Provides additional property $(D flow_call).
template FZbaseInfo(R) {
@ -251,8 +361,7 @@ template FZbaseInfo(R) {
result._flow_call = _current_flow_call;
debug {
if ((_rev && result.base != result._flow_call.base.complement)
|| (!_rev && result.base != result._flow_call.base)) {
if (result.base != result._flow_call.base) {
import std.stdio;
stderr.writeln("invalid flow call at ", _read_name, ": ", result.position);
}
@ -302,28 +411,33 @@ template CIGARbaseInfo(R) {
}
private {
CigarOperation _cigar_operation;
ulong _reference_position;
CigarOperation _cigar_operation = void;
ulong _reference_position = void;
}
}
mixin template rangeMethods() {
private {
const(CigarOperation)[] _cigar = void;
ReversableRange!(identity, const(CigarOperation)[]) _cigar = void;
long _index = void;
CigarOperation _current_cigar_op = void;
ulong _at = void;
ulong _ref_pos = void;
}
/// Current CIGAR operation, available to all extensions
const(CigarOperation) current_cigar_operation() @property const {
return _current_cigar_op;
}
void setup(Args...)(const ref R read, Args)
{
_cigar = read.cigar;
_cigar = reversableRange(read.cigar, read.is_reverse_strand);
_index = _rev ? _cigar.length : -1;
_ref_pos = _rev ? (read.position + read.basesCovered() - 1)
: read.position;
_index = -1;
_ref_pos = reverse_strand ? (read.position + read.basesCovered() - 1)
: read.position;
_moveToNextCigarOperator();
}
@ -337,7 +451,7 @@ template CIGARbaseInfo(R) {
{
++_at;
if (_current_cigar_op.is_reference_consuming) {
_ref_pos += _rev ? -1 : 1;
_ref_pos += reverse_strand ? -1 : 1;
}
if (_at == _current_cigar_op.length) {
@ -355,22 +469,17 @@ template CIGARbaseInfo(R) {
private void _moveToNextCigarOperator() {
_at = 0;
if (!_rev) {
for (++_index; _index < _cigar.length; ++_index)
for (++_index; _index < _cigar.length; ++_index)
{
_current_cigar_op = _cigar[_index];
if (_current_cigar_op.is_query_consuming)
break;
if (_current_cigar_op.is_reference_consuming)
{
_current_cigar_op = _cigar[_index];
if (_current_cigar_op.is_query_consuming)
break;
if (_current_cigar_op.is_reference_consuming)
_ref_pos += _current_cigar_op.length;
}
} else {
for (--_index; _index >= 0; --_index) {
_current_cigar_op = _cigar[_index];
if (_current_cigar_op.is_query_consuming)
break;
if (_current_cigar_op.is_reference_consuming)
if (reverse_strand)
_ref_pos -= _current_cigar_op.length;
else
_ref_pos += _current_cigar_op.length;
}
}
}

115
bio/bam/md/operation.d

@ -1,6 +1,12 @@
module bio.bam.md.operation;
import bio.core.base;
import bio.core.sequence;
import std.conv;
import std.traits;
import std.bitmanip;
import std.algorithm;
/// MD tag operation types
enum MdOperationType : ubyte {
@ -11,35 +17,99 @@ enum MdOperationType : ubyte {
/// Single MD operation.
struct MdOperation {
MdOperationType type; /// Operation type
union {
uint match; /// If this is a match, contains the number of matched bases.
string deletion; /// If this is deletion, contains deleted sequence.
char mismatch; /// If this is a mismatch, contains the mismatched reference base.
alias ReturnType!(nucleotideSequence!SliceableString) DeletionSequence;
private {
MdOperationType _type;
union {
uint _match;
DeletionSequence _deletion;
Base16 _mismatch;
}
}
/// Operation type
MdOperationType type() @property const {
return _type;
}
/// ditto
void type(MdOperationType t) @property {
_type = t;
}
/// Convenience methods
bool is_deletion() @property const {
return _type == MdOperationType.Deletion;
}
/// ditto
bool is_match() @property const {
return _type == MdOperationType.Match;
}
/// ditto
bool is_mismatch() @property const {
return _type == MdOperationType.Mismatch;
}
/// The number of matched bases
ref uint match() @property {
return _match;
}
/// Mismatched reference base
Base16 mismatch() @property const {
return _mismatch;
}
/// ditto
void mismatch(Base16 base) @property {
_mismatch = base;
}
/// Deleted sequence
ref DeletionSequence deletion() @property {
return _deletion;
}
static MdOperation createMatch(uint match) {
MdOperation m;
m.type = MdOperationType.Match;
m.match = match;
m._type = MdOperationType.Match;
m._match = match;
return m;
}
static MdOperation createDeletion(string deletion) {
MdOperation m;
m.type = MdOperationType.Deletion;
m.deletion = deletion;
m._type = MdOperationType.Deletion;
m._deletion = nucleotideSequence(sliceableString(deletion));
return m;
}
static MdOperation createMismatch(char mismatch) {
MdOperation m;
m.type = MdOperationType.Mismatch;
m.mismatch = mismatch;
m._type = MdOperationType.Mismatch;
m._mismatch = Base16(mismatch);
return m;
}
bool opEquals(ref const MdOperation other) const {
static MdOperation createDeletion(DeletionSequence seq) {
MdOperation m;
m._type = MdOperationType.Deletion;
m._deletion = seq;
return m;
}
static MdOperation createMismatch(Base16 base) {
MdOperation m;
m._type = MdOperationType.Mismatch;
m._mismatch = base;
return m;
}
bool opEquals(ref MdOperation other) {
if (type != other.type) {
return false;
@ -47,22 +117,33 @@ struct MdOperation {
final switch (type) {
case MdOperationType.Match:
return match == other.match;
return _match == other._match;
case MdOperationType.Mismatch:
return mismatch == other.mismatch;
case MdOperationType.Deletion:
return deletion == other.deletion;
return equal(_deletion, other._deletion);
}
}
string toString() const {
final switch (type) {
case MdOperationType.Match:
return "Match(" ~ to!string(match) ~ ")";
return "Match(" ~ to!string(_match) ~ ")";
case MdOperationType.Mismatch:
return "Mismatch(" ~ to!string(mismatch) ~ ")";
return "Mismatch(" ~ to!string(_mismatch) ~ ")";
case MdOperationType.Deletion:
return "Deletion(" ~ to!string(deletion) ~ ")";
return "Deletion(" ~ to!string(_deletion) ~ ")";
}
}
}
/// Returns MD operation with reverse-complemented data
MdOperation reverseMdOp(MdOperation op) {
if (op.is_deletion)
return MdOperation.createDeletion(op.deletion.reverse);
if (op.is_mismatch)
return MdOperation.createMismatch(op.mismatch.complement);
return op;
}

139
bio/bam/md/parse.d

@ -5,56 +5,132 @@ import std.ascii;
import std.array;
import std.algorithm;
import std.functional;
import std.range;
import std.conv;
import std.traits;
/// Returns range of MD operations. Zero matches are skipped.
/// Returns bidirectional range of MD operations. Zero matches are skipped.
auto mdOperations(string md) {
struct Result {
private {
MdOperation _cur_md_op; // current MD operation
bool _empty = false;
string _md;
string _md = void;
MdOperation _cached_front = void;
MdOperation _cached_back = void;
ubyte _rem = 255;
}
this(string md) {
_md = md;
popFront();
if (!cacheFront()) {
_rem = 0;
} else {
if (!cacheBack()) {
_cached_back = _cached_front;
_rem = 1;
}
}
}
bool empty() @property {
return _empty;
return _rem == 0;
}
Result save() @property {
Result res = void;
res._md = _md;
res._cached_front = _cached_front;
res._cached_back = _cached_back;
res._rem = _rem;
return res;
}
ref MdOperation front() @property {
return _cached_front;
}
MdOperation front() @property {
return _cur_md_op;
ref MdOperation back() @property {
return _cached_back;
}
void popFront() {
if (_md.empty) {
_empty = true;
return;
if (_rem == 255) {
_cached_front = _cached_back;
_rem = 1;
} else {
_rem = 0;
}
} else {
if (!cacheFront())
_rem = 0;
}
}
if (_md[0] == '^') {
void popBack() {
if (_md.empty) {
if (_rem == 255) {
_cached_back = _cached_front;
_rem = 1;
} else {
_rem = 0;
}
} else {
if (!cacheBack())
_rem = 0;
}
}
private bool cacheFront() {
if (_md.empty)
return false;
if (_md[0] == '^') { // deletion, get bases
_md = _md[1 .. $];
auto len = countUntil!(not!isUpper)(_md);
if (len == -1) {
len = _md.length;
}
_cur_md_op = MdOperation.createDeletion(_md[0 .. len]);
_cached_front = MdOperation.createDeletion(_md[0 .. len]);
_md = _md[len .. $];
} else if (isDigit(_md[0])) {
} else if (isDigit(_md[0])) { // match, get number
auto len = countUntil!(not!isDigit)(_md);
if (len == -1) {
len = _md.length;
}
_cur_md_op = MdOperation.createMatch(to!uint(_md[0 .. len]));
_cached_front = MdOperation.createMatch(to!uint(_md[0 .. len]));
_md = _md[len .. $];
} else {
_cur_md_op = MdOperation.createMismatch(_md[0]);
} else { // mismatch
_cached_front = MdOperation.createMismatch(_md[0]);
_md = _md[1 .. $];
}
return true;
}
private bool cacheBack() {
if (_md.empty)
return false;
if (isDigit(_md[$ - 1])) { // match, get number
auto len = countUntil!(not!isDigit)(retro(_md));
if (len == -1) {
len = _md.length;
}
_cached_back = MdOperation.createMatch(to!uint(_md[$ - len .. $]));
_md = _md[0 .. $ - len];
} else {
if (_md.length == 1 || isDigit(_md[$ - 2])) { // mismatch
_cached_back = MdOperation.createMismatch(_md[$ - 1]);
_md = _md[0 .. $ - 1];
} else { // deletion
auto len = countUntil!"a == '^'"(retro(_md));
_cached_back = MdOperation.createDeletion(_md[$ - len .. $]);
_md = _md[0 .. $ - len - 1];
}
}
return true;
}
}
@ -63,13 +139,18 @@ auto mdOperations(string md) {
op.match == 0;
}
return filter!(not!isZeroMatch)(Result(md));
return filterBidirectional!(not!isZeroMatch)(Result(md));
}
/// Alias for return type of mdOperations
alias ReturnType!mdOperations MdOperationRange;
unittest {
import std.algorithm;
import std.stdio;
assert(equal(mdOperations("86"),
[MdOperation.createMatch(86)]));
@ -82,17 +163,17 @@ unittest {
MdOperation.createDeletion("T"),
MdOperation.createMatch(28)]));
assert(equal(mdOperations("3C6C0A13^A4C2"),
[MdOperation.createMatch(3),
MdOperation.createMismatch('C'),
MdOperation.createMatch(6),
MdOperation.createMismatch('C'),
MdOperation.createMismatch('A'),
MdOperation.createMatch(13),
MdOperation.createDeletion("A"),
MdOperation.createMatch(4),
MdOperation.createMismatch('C'),
MdOperation.createMatch(2)]));
assert(equal(retro(mdOperations("3C6C0A13^A4C2")),
retro([MdOperation.createMatch(3),
MdOperation.createMismatch('C'),
MdOperation.createMatch(6),
MdOperation.createMismatch('C'),
MdOperation.createMismatch('A'),
MdOperation.createMatch(13),
MdOperation.createDeletion("A"),
MdOperation.createMatch(4),
MdOperation.createMismatch('C'),
MdOperation.createMatch(2)])));
assert(equal(mdOperations("27^TTT63"),
[MdOperation.createMatch(27),

2
bio/bam/md/reconstruct.d

@ -159,7 +159,7 @@ auto dna(T)(T read)
}
break;
case MdOperationType.Deletion:
_cur_md_op.deletion = _cur_md_op.deletion[1..$];
_cur_md_op.deletion.popFront();
if (_cur_md_op.deletion.empty) {
_fetchNextMdOperation();
}

30
bio/core/base.d

@ -231,3 +231,33 @@ unittest {
assert((cast(Base16)b5).internal_code == 8);
}
/// Complement base, which might be Base5, Base16, char, or dchar.
B complementBase(B)(B base) {
static if(is(Unqual!B == dchar) || is(Unqual!B == char))
{
return cast(B)(Base16(base).complement);
}
else
return base.complement;
}
/// Convert character to base
template charToBase(B=Base16)
{
B charToBase(C)(C c)
if(is(Unqual!C == char) || is(Unqual!C == dchar))
{
return B(c);
}
}
unittest {
assert(complementBase('T') == 'A');
assert(complementBase('G') == 'C');
assert(complementBase(Base5('A')) == Base5('T'));
assert(complementBase(Base16('C')) == Base16('G'));
assert(charToBase!Base16('A').complement == Base16('T'));
}

223
bio/core/sequence.d

@ -0,0 +1,223 @@
/*
This file is part of BioD.
Copyright (C) 2012 Artem Tarasov <lomereiter@gmail.com>
BioD is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
BioD is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
module bio.core.sequence;
import bio.core.base;
import std.algorithm;
import std.range;
import std.conv;
import std.traits;
public import std.array;
/// Identity function
T identity(T)(auto ref T t) { return t; }
/// Range that allows to unify operations in forward and reverse directions
/// without virtual function call overhead introduced by $(D inputRangeObject).
///
/// $(D reverseTransform) is a function that will be applied to elements
/// if range is iterated backwards.
struct ReversableRange(alias reverseTransform=identity, R)
if(isBidirectionalRange!R)
{
private
{
R _range = void;
bool _rev = void;
}
/// Construct reversable range.
///
/// Params:
/// range = bidirectional range
/// reverse = if true, all operations on the range will be as if
/// $(D retro(range)) was used instead of $(D range).
this(R range, bool reverse=false)
{
_rev = reverse;
_range = range;
}
/// Bidirectional range primitives
bool empty() @property
{
return _range.empty;
}
/// ditto
auto front() @property
{
return _rev ? reverseTransform(_range.back) : _range.front;
}
/// ditto
auto back() @property
{
return _rev ? reverseTransform(_range.front) : _range.back;
}
/// ditto
void popFront()
{
if (_rev)
_range.popBack();
else
_range.popFront();
}
/// ditto
void popBack()
{
if (_rev)
_range.popFront();
else
_range.popBack();
}
/// ditto
auto save() @property
{
return ReversableRange(_range.save, _rev);
}
/// Reverse of this range
ReversableRange reverse() @property {
return ReversableRange(_range.save, !_rev);
}
static if(hasLength!R)
{
/// If source range has length, the result also has length
size_t length() @property
{
return _range.length;
}
}
static if(isRandomAccessRange!R)
{
/// If source range is a random access range, $(D opIndex) is defined
auto opIndex(size_t index)
{
if (_rev)
return reverseTransform(_range[_range.length - 1 - index]);
else
return _range[index];
}
}
static if(hasSlicing!R)
{
/// Slicing is also propagated
auto opSlice(size_t from, size_t to)
{
if (_rev)
{
auto len = _range.length;
//
// [b, e) -> (len - 1 - e, len - 1 - b] ~ [len - e, len - b)
//
return ReversableRange(_range[len - to .. len - from], true);
}
else
return ReversableRange(_range[from .. to], false);
}
}
}
/// Create reversable range from bidirectional one.
ReversableRange!(reverseTransform, R)
reversableRange(alias reverseTransform=identity, R)(R range, bool reverse=false)
{
return typeof(return)(range, reverse);
}
unittest {
auto bidir_range = [1, 2, 3, 4, 5];
auto rev = reversableRange(bidir_range[], true);
assert(rev.front == 5);
assert(rev[2] == 3);
rev.popFront();
assert(rev.back == 1);
assert(rev.front == 4);
assert(equal(rev[1 .. 3], [3, 2]));
// Here. That's the whole point.
// One can't do the same with $(D retro)
// without using $(D inputRangeObject),
// but that kills performance because
// virtual calls can not be inlined.
rev = reversableRange(bidir_range[], false);
assert(rev.front == 1);
assert(equal(rev[1 .. 3], [2, 3]));
}
/// Sequence of bases. Element of reversed range will be complemented.
template Sequence(R)
{
alias ReversableRange!(complementBase, R) Sequence;
}
/// Returns an object very similar to string, but sliceable.
/// Tricks std.traits.isNarrowString.
auto sliceableString(string s) {
return map!"cast(char)a"(cast(ubyte[])s);
}
///
alias ReturnType!sliceableString SliceableString;
/// Create nucleotide sequence from bidirectional base range.
auto nucleotideSequence(R)(R bases, bool reverse=false)
if(isBidirectionalRange!R)
{
static if(isNarrowString!R)
{
return nucleotideSequence(sliceableString(bases), reverse);
}
else static if(is(Unqual!(ElementType!R) == char) ||
is(Unqual!(ElementType!R) == dchar))
{
return nucleotideSequence(map!(charToBase!Base16)(bases), reverse);
}
else
{
return Sequence!R(bases, reverse);
}
}
unittest {
auto seq0 = nucleotideSequence("ACGTACGT");
// reverse-complement
assert(equal(seq0.reverse[2 .. 6], "GTAC"));
auto seq1 = nucleotideSequence(seq0, true);
assert(equal(seq1[1 .. 5], "CGTA"));
assert(equal(seq1, map!complementBase(retro(seq0))));
seq1 = nucleotideSequence(seq0, false);
assert(equal(seq1, seq0));
}

24
test/unittests.d

@ -36,6 +36,10 @@ import bio.core.utils.roundbuf;
import bio.core.utils.range;
import bio.core.utils.tmpfile;
import bio.core.utils.stream;
import bio.core.sequence;
import bio.core.base;
import bio.core.tinymap;
import bio.core.utils.roundbuf;
import std.path;
import std.range;
@ -350,10 +354,14 @@ unittest {
auto read = reads[1];
assert(!read.is_reverse_strand);
auto basesFZ = basesWith!"FZ"(read, arg!"flowOrder"(flow_order),
auto bases = basesWith!("FZ", "MD")(read,
arg!"flowOrder"(flow_order),
arg!"keySequence"(key_sequence));
assert(equal(basesFZ.save, read.sequence));
assert(equal(take(map!"a.flow_call.intensity_value"(basesFZ.save), 92),
assert(equal(drop(map!"a.reference_base"(bases.save), 191),
"-CCCGATTGGTCGTTGCTTTACGCTGATTGGCGAGTCCGGGGAACGTACCTTTGCTATCAGTCCAGGCCACATGAACCAGCTGCGGGCTGAAAGCATTCCGGAAGATGTGATTGCCGGACCTCGGCACTGGTTCTCACCTCATATCTGGTGCGTTGCAAGCCGGGTGAACCCATGCCGGAAGCACCATGAAAGCCATTGAGTACGCGAAGAAATATA"));
assert(equal(bases.save, read.sequence));
assert(equal(take(map!"a.flow_call.intensity_value"(bases.save), 92),
[219, 219, 194, 194, 92, 107, 83, 198, 198, 78,
// A A C C T G A T T A
292, 292, 292, 81, 79, 78, 95, 99, 315, 315, 315,
@ -377,12 +385,14 @@ unittest {
if (r.is_unmapped) continue;
if (r.cigar.length == 0) continue;
if (r.is_reverse_strand) {
basesFZ = basesWith!"FZ"(r, arg!"flowOrder"(flow_order), arg!"keySequence"(key_sequence));
bases = basesWith!("FZ", "MD")(r, arg!"flowOrder"(flow_order),
arg!"keySequence"(key_sequence));
// if reverse strand, bases are also reverse complemented
assert(equal(basesFZ.save, map!"a.complement"(retro(r.sequence))));
assert(equal(bases.save, map!"a.complement"(retro(r.sequence))));
} else {
basesFZ = basesWith!"FZ"(r, arg!"flowOrder"(flow_order), arg!"keySequence"(key_sequence));
assert(equal(basesFZ.save, r.sequence));
bases = basesWith!("FZ", "MD")(r, arg!"flowOrder"(flow_order),
arg!"keySequence"(key_sequence));
assert(equal(bases.save, r.sequence));
}
}
}

Loading…
Cancel
Save