lib.py 66 KB
Newer Older
jhoogenboom's avatar
jhoogenboom committed
1 2
#!/usr/bin/env python

3
import re, sys, argparse, random, itertools
jhoogenboom's avatar
jhoogenboom committed
4
#import numpy as np  # Imported only when calling nnls()
jhoogenboom's avatar
jhoogenboom committed
5 6

from ConfigParser import RawConfigParser, MissingSectionHeaderError
7
from StringIO import StringIO
jhoogenboom's avatar
jhoogenboom committed
8 9 10 11

# Patterns that match entire sequences.
PAT_SEQ_RAW = re.compile("^[ACGT]*$")
PAT_SEQ_TSSV = re.compile("^(?:[ACGT]+\(\d+\))*$")
12
PAT_SEQ_ALLELENAME_STR = re.compile(  # First line: n_ACT[m] or alias.
jhoogenboom's avatar
jhoogenboom committed
13
    "^(?:(?:(?:CE)?-?\d+(?:\.\d+)?_(?:[ACGT]+\[\d+\])*)|((?!_).+?))"
jhoogenboom's avatar
jhoogenboom committed
14
    "(?:_[-+]\d+(?:\.1)?(?P<a>(?:(?<=\.1)-)|(?<!\.1)[ACGT]+)>"  # _+3A>
15
        "(?!(?P=a))(?:[ACGT]+|-))*$")  # Portion of variants after '>'.
16 17 18
PAT_SEQ_ALLELENAME_SNP = re.compile(
    "^REF$|^(?:(?:(?<=^)|(?<!^) )"  # 'REF' or space-separated variants.
    "\d+(?:\.1)?(?P<a>(?:(?<=\.1)-)|(?<!\.1)[ACGT]+)>"
19
        "(?!(?P=a))(?:[ACGT]+|-))+$")  # Portion of variants after '>'.
20 21
PAT_SEQ_ALLELENAME_MT = re.compile(
    "^REF$|^(?:(?:(?<=^)|(?<!^) )"  # 'REF' or space-separated variants.
22
    "(?:-?\d+\.\d+[ACGT]|(?P<a>[ACGT])?\d+(?(a)(?!(?P=a)))(?:[ACGT-]|del)))+$")
jhoogenboom's avatar
jhoogenboom committed
23

24
# Patterns that match blocks of TSSV-style sequences and allele names.
jhoogenboom's avatar
jhoogenboom committed
25 26
PAT_TSSV_BLOCK = re.compile("([ACGT]+)\((\d+)\)")
PAT_ALLELENAME_BLOCK = re.compile("([ACGT]+)\[(\d+)\]")
27
PAT_ALIAS = re.compile("^(?!_).+$")
jhoogenboom's avatar
jhoogenboom committed
28

29 30 31 32
# Patterns that match a single variant.
PAT_VARIANT_STR = re.compile(
    "^(?P<pos>[-+]\d+)(?:\.(?P<ins>1))?"
    "(?P<old>(?:(?<=\.1)-)|(?<!\.1)[ACGT]+)>"
33
    "(?!(?P=old))(?P<new>[ACGT]+|-)$")
34 35 36
PAT_VARIANT_SNP = re.compile(
    "^(?P<pos>\d+)(?:\.(?P<ins>1))?"
    "(?P<old>(?:(?<=\.1)-)|(?<!\.1)[ACGT]+)>"
37
    "(?!(?P=old))(?P<new>[ACGT]+|-)$")
38 39 40 41
PAT_VARIANT_MT = re.compile(
    "^(?P<old>(?P<a>[ACGT])|-?)"
    "(?P<pos>\d+)(?(a)|(?:\.(?P<ins>\d+))?)"
    "(?P<new>[ACGT-]|del)$")
jhoogenboom's avatar
jhoogenboom committed
42 43

# Patterns that match (parts of) an STR definition.
44
PAT_STR_DEF = re.compile("^(?:(?:(?<=^)|(?<!^)\s+)[ACGT]+\s+\d+\s+\d+)*$")
jhoogenboom's avatar
jhoogenboom committed
45 46 47
PAT_STR_DEF_BLOCK = re.compile("([ACGT]+)\s+(\d+)\s+(\d+)")

# Pattern to split a comma-, semicolon-, or space-separated list.
48 49 50 51 52
PAT_SPLIT = re.compile("\s*[,; \t]\s*")

# Pattern that matches a chromosome name/number.
PAT_CHROMOSOME = re.compile(
    "^(?:[Cc][Hh][Rr](?:[Oo][Mm])?)?([1-9XYM]|1\d|2[0-2])$")
jhoogenboom's avatar
jhoogenboom committed
53

jhoogenboom's avatar
jhoogenboom committed
54 55
# Default regular expression to capture sample tags in file names.
# This is the default of the -e command line option.
jhoogenboom's avatar
jhoogenboom committed
56
DEF_TAG_EXPR = "^(.*?)(?:\.[^.]+)?$"
jhoogenboom's avatar
jhoogenboom committed
57 58 59 60 61

# Default formatting template to write sample tags.
# This is the default of the -f command line option.
DEF_TAG_FORMAT = "\\1"

62 63 64 65 66
# Default formatting template to construct output file names for batch
# processing.  \1 and \2 refer to sample tag and tool name.
# This is the default for the -o command line option with batch support.
DEF_OUTFILE_FORMAT = "\\1-\\2.out"

jhoogenboom's avatar
jhoogenboom committed
67 68 69 70 71 72
# IUPAC Table of complementary bases.
COMPL = {"A": "T", "T": "A", "U": "A", "G": "C", "C": "G", "R": "Y", "Y": "R",
         "K": "M", "M": "K", "B": "V", "V": "B", "D": "H", "H": "D",
         "a": "t", "t": "a", "u": "a", "g": "c", "c": "g", "r": "y", "y": "r",
         "k": "m", "m": "k", "b": "v", "v": "b", "d": "h", "h": "d"}

73 74 75
# Special values that may appear in the place of a sequence.
SEQ_SPECIAL_VALUES = ("No data", "Other sequences")

jhoogenboom's avatar
jhoogenboom committed
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
def get_genome_pos(location, x, invert=False):
    """Get the genome position of the x-th base in a sequence."""
    if invert:
        offset = 0
        for i in range(1, len(location)):
            if i % 2:
                # Starting position.
                pos = location[i]
            elif pos <= x <= location[i]:
                # x is in the current range
                break
            else:
                offset += location[i]-pos+1
        else:
            if len(location) % 2:
                raise ValueError("Position %i is outside sequence range" % x)
        return offset + x - pos
    else:
        for i in range(1, len(location)):
            if i % 2:
                # Starting position.
                pos = location[i]
            elif location[i]-pos < x:
                # x is after this ending position
                x -= location[i]-pos+1
            else:
                # x is before this ending position
                break
        return pos + x
#get_genome_pos


109
def call_variants(template, sequence, location="suffix", cache=True,
jhoogenboom's avatar
jhoogenboom committed
110 111 112
                  debug=False):
    """
    Perform a global alignment of sequence to template and return a
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
    list of variants detected.  The format (nomenclature) of the
    returned variants depends on the location argument.

    If location is "suffix" (the default), all variants are given as
    substitutions in the form posX>Y, where the first base in the
    template is pos=1.  With location set to "prefix", bases are counted
    from right to left instead.  Insertions and deletions are written as
    pos.1->Y and posX>-, respectively.

    If location is a tuple ("M", position) with any integer for the
    position, variants are written following the mtDNA nomenclature
    guidelines.  The given position is that of the first base in the
    template.

    If location is a tuple ("chromosome name", position), a
    NotImplementedError is raised.
jhoogenboom's avatar
jhoogenboom committed
129 130 131 132 133 134 135 136 137

    By default, the results of this function are cached.  Set cache to
    False to suppress caching the result and reduce memory usage.

    Setting debug to True will cause the alignment matrices to be
    printed to sys.stdout.  Be aware that they can be quite large.
    """
    # Saving the results in a cache to avoid repeating alignments.
    try:
138
        return call_variants.cache[template, sequence, location]
jhoogenboom's avatar
jhoogenboom committed
139
    except KeyError:
140
        cache_key = location
jhoogenboom's avatar
jhoogenboom committed
141 142 143 144 145

    row_offset = len(template) + 1
    matrix_match = [0] * row_offset * (len(sequence)+1)
    matrix_gap1 = [-sys.maxint-1] * row_offset * (len(sequence)+1)
    matrix_gap2 = [-sys.maxint-1] * row_offset * (len(sequence)+1)
146 147 148 149 150 151 152 153 154 155 156 157 158
    matrix_direction = [0] * row_offset * (len(sequence)+1)

    # Matrix and arrow enum constants.
    M_MATCH = 0
    M_GAP1 = 1
    M_GAP2 = 2
    A_MATCH = 1
    A_HORZ_O = 2
    A_HORZ_E = 4
    A_VERT_O = 8
    A_VERT_E = 16

    # Settings.
jhoogenboom's avatar
jhoogenboom committed
159
    MATCH_SCORE = 1
160 161 162
    MISMATCH_SCORE = -3
    GAP_OPEN_SCORE = -7
    GAP_EXTEND_SCORE = -2
163 164 165 166 167 168 169 170
    variant_format = "%i%s>%s"

    if location == "prefix":
        location = ("prefix", -len(template))
    elif location == "suffix":
        # Include plus signs for position numbers.
        variant_format = "%+i%s>%s"
        location = ("suffix", 1)
171
    elif type(location) != tuple or len(location) < 2:
172
        raise ValueError("Unknown location %r. It should be 'prefix', "
173 174
            "'suffix', or a tuple (chromosome, position [, endpos])" %
            location)
jhoogenboom's avatar
jhoogenboom committed
175 176 177 178 179 180 181 182 183 184 185

    for i in range(len(matrix_match)):
        x = i % row_offset
        y = i / row_offset

        # Initialisation of first row and column.
        if x == 0 or y == 0:
            if x != 0:
                # Top row.
                matrix_gap1[i] = GAP_OPEN_SCORE + GAP_EXTEND_SCORE * (x - 1)
                matrix_match[i] = matrix_gap1[i]
186
                matrix_direction[i] = A_HORZ_E | (A_HORZ_O if x == 1 else 0)
jhoogenboom's avatar
jhoogenboom committed
187 188 189 190
            elif y != 0:
                # Left column.
                matrix_gap2[i] = GAP_OPEN_SCORE + GAP_EXTEND_SCORE * (y - 1)
                matrix_match[i] = matrix_gap2[i]
191 192 193 194
                matrix_direction[i] = A_VERT_E | (A_VERT_O if y == 1 else 0)
            else:
                # Top left corner.
                matrix_direction[i] = A_MATCH
jhoogenboom's avatar
jhoogenboom committed
195 196 197 198 199 200 201
            continue

        if template[x-1] == sequence[y-1]:
            match = MATCH_SCORE
        else:
            match = MISMATCH_SCORE

202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
        options_gap1 = (
            matrix_match[i-1] + GAP_OPEN_SCORE,
            matrix_gap1[i-1] + GAP_EXTEND_SCORE)
        matrix_gap1[i] = max(options_gap1)
        if options_gap1[0] > options_gap1[1]:
            matrix_direction[i] |= A_HORZ_O  # Must exit M_GAP1 here.

        options_gap2 = (
            matrix_match[i-row_offset] + GAP_OPEN_SCORE,
            matrix_gap2[i-row_offset] + GAP_EXTEND_SCORE)
        matrix_gap2[i] = max(options_gap2)
        if options_gap2[0] > options_gap2[1]:
            matrix_direction[i] |= A_VERT_O  # Must exit M_GAP2 here.

        options = (
jhoogenboom's avatar
jhoogenboom committed
217 218 219
            matrix_match[i-1-row_offset] + match,
            matrix_gap1[i],
            matrix_gap2[i])
220 221 222 223 224 225 226
        matrix_match[i] = max(options)
        if options[0] == matrix_match[i]:
            matrix_direction[i] |= A_MATCH  # Can stay in M_MATCH here.
        if options[1] == matrix_match[i]:
            matrix_direction[i] |= A_HORZ_E  # Can enter M_GAP1 here.
        if options[2] == matrix_match[i]:
            matrix_direction[i] |= A_VERT_E  # Can enter M_GAP2 here.
jhoogenboom's avatar
jhoogenboom committed
227 228 229 230 231 232 233 234 235 236 237

    if debug:
        print("GAP1")
        for i in range(0, len(matrix_gap1), row_offset):
            print(("%5i" * row_offset) % tuple(matrix_gap1[i:i+row_offset]))
        print("GAP2")
        for i in range(0, len(matrix_gap2), row_offset):
            print(("%5i" * row_offset) % tuple(matrix_gap2[i:i+row_offset]))
        print("Match")
        for i in range(0, len(matrix_match), row_offset):
            print(("%5i" * row_offset) % tuple(matrix_match[i:i+row_offset]))
238 239 240 241 242 243 244 245 246 247
        print("FLAGS")
        for i in range(0, len(matrix_direction), row_offset):
            print(("%5s|" * row_offset) % tuple("".join([
                "h" if x & A_HORZ_O else " ",
                "H" if x & A_HORZ_E else " ",
                "D" if x & A_MATCH else " ",
                "V" if x & A_VERT_E else " ",
                "v" if x & A_VERT_O else " "
            ]) for x in matrix_direction[i:i+row_offset]))
        print("Traceback")
jhoogenboom's avatar
jhoogenboom committed
248 249 250 251 252 253 254


    # Backtracking.
    variants = []
    variant_template = 0
    variant_sequence = 0
    i = len(matrix_match) - 1
255
    in_matrix = M_MATCH  # May change before first step.
jhoogenboom's avatar
jhoogenboom committed
256 257 258
    while i >= 0:
        x = i % row_offset
        y = i / row_offset
259 260 261 262 263 264 265 266 267 268 269 270 271 272
        if debug:
            print("(%i, %i)" % (x,y))

        if in_matrix == M_MATCH:
            # Make gaps as soon as possible (pushed right).
            if matrix_direction[i] & A_HORZ_E:
                in_matrix = M_GAP1
            elif matrix_direction[i] & A_VERT_E:
                in_matrix = M_GAP2
            elif not (matrix_direction[i] & A_MATCH):
                raise ValueError(
                    "Alignment error: Dead route! (This is a bug.) [%s,%s]" % (template,sequence))

        if in_matrix == M_GAP1:
jhoogenboom's avatar
jhoogenboom committed
273 274
            # Go horizontally.  Deletion.
            variant_template += 1
275 276 277
            if matrix_direction[i] & A_HORZ_O:
                # End of gap, go diagonally after this.
                in_matrix = M_MATCH
jhoogenboom's avatar
jhoogenboom committed
278 279 280
            i -= 1
            continue

281
        if in_matrix == M_GAP2:
jhoogenboom's avatar
jhoogenboom committed
282 283
            # Go vertically.  Insertion.
            variant_sequence += 1
284 285 286
            if matrix_direction[i] & A_VERT_O:
                # End of gap, go diagonally after this.
                in_matrix = M_MATCH
jhoogenboom's avatar
jhoogenboom committed
287 288 289 290
            i -= row_offset
            continue

        # Go diagonally.  Either match or mismatch.
291 292 293 294 295 296
        if i != 0 and template[x - 1] != sequence[y - 1]:
            # Start/extend mismatch.
            variant_template += 1
            variant_sequence += 1

        else:
jhoogenboom's avatar
jhoogenboom committed
297 298
            # Match.  Flush variants.
            if variant_template or variant_sequence:
299 300 301 302 303 304
                if location[0] == "M":
                    # MtDNA variants are one-base-at-a-time.
                    for j in range(
                            max(variant_template, variant_sequence)-1, -1, -1):
                        variants.append("%s%i%s%s" % (
                            template[x+j] if j < variant_template else "",#"-",
305 306
                            get_genome_pos(
                                location, x + min(j, variant_template-1)),
307 308
                            ".%i" % (j-variant_template+1)
                                if j >= variant_template else "",
309
                            sequence[y+j] if j < variant_sequence else "del"))
310
                elif variant_template == 0:
jhoogenboom's avatar
jhoogenboom committed
311
                    # Insertions: "-131.1->C" instead of "-130->C".
312
                    variants.append(variant_format % (
313
                        get_genome_pos(location, x - 1),
314
                        ".1-",
jhoogenboom's avatar
jhoogenboom committed
315 316
                        sequence[y:y+variant_sequence]))
                else:
317
                    variants.append(variant_format % (
318
                        get_genome_pos(location, x),
jhoogenboom's avatar
jhoogenboom committed
319 320 321 322 323 324
                        template[x:x+variant_template],
                        sequence[y:y+variant_sequence] or "-"))
                variant_template = 0
                variant_sequence = 0
        i -= 1 + row_offset

325 326
    # Variants were called from right to left.  Reverse their order.
    if location[0] != "prefix":
jhoogenboom's avatar
jhoogenboom committed
327 328 329 330
        variants.reverse()

    # Store the result in the cache.
    if cache:
331
        call_variants.cache[template, sequence, cache_key] = variants
jhoogenboom's avatar
jhoogenboom committed
332 333 334 335 336
    return variants
#call_variants
call_variants.cache = {}


337
def mutate_sequence(seq, variants, location=None):
jhoogenboom's avatar
jhoogenboom committed
338
    """Apply the given variants to the given sequence."""
339
    if type(location) != tuple or len(location) < 2:
340
        pattern = PAT_VARIANT_STR
341
        location = (None, 0)
342 343
    elif location[0] == "M":
        pattern = PAT_VARIANT_MT
344
        location = (location[0], location[1]-1) + tuple(location[2:])
345 346
    else:
        pattern = PAT_VARIANT_SNP
347
        location = (location[0], location[1]-1) + tuple(location[2:])
348 349

    seq = [[]] + [[base] for base in seq]
jhoogenboom's avatar
jhoogenboom committed
350
    for variant in variants:
351
        vm = pattern.match(variant)
jhoogenboom's avatar
jhoogenboom committed
352 353
        if vm is None:
            raise ValueError("Unrecognised variant '%s'" % variant)
354 355 356 357
        pos = int(vm.group("pos"))
        ins = int(vm.group("ins") or 0)
        old = vm.group("old")
        new = vm.group("new")
jhoogenboom's avatar
jhoogenboom committed
358 359
        if old == "-":
            old = ""
360
        if new == "-" or new == "del":
jhoogenboom's avatar
jhoogenboom committed
361 362
            new = ""
        if pos < 0:
363
            pos += len(seq)
364 365
        pos = get_genome_pos(location, pos, True)
        if pos < 0 or (pos == 0 and not ins) or pos >= len(seq):
366
            raise ValueError(
367 368
                "Position of variant '%s' is outside sequence range" %
                    (variant))
369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
        if (not ins and old and old != "".join("".join(x[:1])
                for x in seq[pos:pos+len(old)])):
            raise ValueError(
                "Incorrect original sequence in variant '%s'; should be '%s'!"
                % (variant, "".join("".join(x[:1])
                    for x in seq[pos:pos+len(old)])))
        elif not ins and not old:
            # MtDNA substitution with reference base omitted.
            old = "".join("".join(x[:1]) for x in seq[pos:pos+len(new)])
        if not ins:
            # Remove old bases, retaining those inserted between/after.
            seq[pos:pos+len(old)] = [
                [""] + x[1:] for x in seq[pos:pos+len(old)]]
            # Place new entirely in the position of the first old base.
            seq[pos][0] = new
jhoogenboom's avatar
jhoogenboom committed
384
        else:
385 386 387 388 389
            # Insert new exactly ins positions after pos.
            while len(seq[pos]) <= ins:
                seq[pos].append("")
            seq[pos][ins] = new
    return "".join("".join(x) for x in seq)
jhoogenboom's avatar
jhoogenboom committed
390 391 392
#mutate_sequence


393
def parse_library(libfile, stream=False):
jhoogenboom's avatar
jhoogenboom committed
394
    try:
395
        if not stream:
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410
            libfile = sys.stdin if libfile == "-" else open(libfile, "r")
        if libfile == sys.stdin:
            # Can't seek on pipes, so read it into a buffer first.
            libfile = StringIO(sys.stdin.read())
        try:
            library = parse_library_ini(libfile)
            if not stream:
                libfile.close()
            return library
        except MissingSectionHeaderError:
            # Not an ini file.
            pass
        libfile.seek(0)
        library = parse_library_tsv(libfile)
        if not stream and libfile != sys.stdin:
411 412
            libfile.close()
        return library
413 414
    except ValueError as err:
        raise argparse.ArgumentTypeError(err)
jhoogenboom's avatar
jhoogenboom committed
415 416 417 418 419 420 421 422 423 424 425 426 427 428
#parse_library


def parse_library_tsv(handle):
    """
    Parse a TSSV library file (tab-separated values format).

    The provided file should contain at least four columns: marker name,
    left flanking sequence, right flanking sequence, and STR definition.

    Return a nested dict with top-level keys "flanks" and "regex".
    """
    library = {
      "flanks": {},
429
      "regex": {},
430
      "blocks_middle": {}
jhoogenboom's avatar
jhoogenboom committed
431 432
    }
    for line in handle:
jhoogenboom's avatar
jhoogenboom committed
433
        line = [x.strip() for x in line.rstrip("\r\n").split("\t")]
jhoogenboom's avatar
jhoogenboom committed
434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450
        if line == [""]:
            continue
        if len(line) < 4:
            raise ValueError(
                "Invalid library file: encountered line with %i columns, "
                "need at least 4" % len(line))
        marker = line[0]
        if PAT_SEQ_RAW.match(line[1]) is None:
            raise ValueError("Flanking sequence '%s' of marker %s is invalid" %
                             (line[1], marker))
        if PAT_SEQ_RAW.match(line[2]) is None:
            raise ValueError("Flanking sequence '%s' of marker %s is invalid" %
                             (line[2], marker))
        if PAT_STR_DEF.match(line[3]) is None:
            raise ValueError("STR definition '%s' of marker %s is invalid" %
                             (line[3], marker))
        library["flanks"][marker] = line[1:3]
451 452 453
        library["blocks_middle"][marker] = [
            (block[0], int(block[1]), int(block[2])) for block in
                PAT_STR_DEF_BLOCK.findall(line[3])]
454
        library["regex"][marker] = re.compile(
455 456 457
            "".join(("^", "".join(
                "(%s){%s,%s}" % x for x in PAT_STR_DEF_BLOCK.findall(line[3])),
                "$")))
jhoogenboom's avatar
jhoogenboom committed
458 459 460 461 462 463 464 465 466 467
    return library
#parse_library_tsv


def parse_library_ini(handle):
    library = {
      "flanks": {},
      "prefix": {},
      "suffix": {},
      "regex": {},
468
      "blocks_middle": {},
469 470
      "nostr_reference": {},
      "genome_position": {},
471 472
      "length_adjust": {},
      "block_length": {},
473
      "max_expected_copies": {},
474
      "expected_length": {},
jhoogenboom's avatar
jhoogenboom committed
475 476 477 478 479 480 481 482 483 484
      "aliases": {}
    }
    markers = set()

    ini = RawConfigParser()
    ini.optionxform = str
    ini.readfp(handle)
    for section in ini.sections():
        for marker in ini.options(section):
            value = ini.get(section, marker)
485 486
            section_low = section.lower()
            if section_low == "flanks":
jhoogenboom's avatar
jhoogenboom committed
487 488 489 490 491 492 493 494 495 496 497 498
                values = PAT_SPLIT.split(value)
                if len(values) != 2:
                    raise ValueError(
                        "For marker %s, %i flanking sequences were given,"
                        "need exactly 2" % (marker, len(values)))
                for value in values:
                    if PAT_SEQ_RAW.match(value) is None:
                        raise ValueError(
                            "Flanking sequence '%s' of marker %s is invalid" %
                            (value, marker))
                library["flanks"][marker] = values
                markers.add(marker)
499
            elif section_low == "prefix":
500 501 502
                if marker in library["nostr_reference"]:
                    raise ValueError(
                        "A prefix was defined for non-STR marker %s" % marker)
jhoogenboom's avatar
jhoogenboom committed
503 504 505 506 507 508
                values = PAT_SPLIT.split(value)
                for value in values:
                    if PAT_SEQ_RAW.match(value) is None:
                        raise ValueError(
                            "Prefix sequence '%s' of marker %s is invalid" %
                            (value, marker))
509
                library["prefix"][marker] = values
jhoogenboom's avatar
jhoogenboom committed
510
                markers.add(marker)
511
            elif section_low == "suffix":
512 513 514
                if marker in library["nostr_reference"]:
                    raise ValueError(
                        "A suffix was defined for non-STR marker %s" % marker)
jhoogenboom's avatar
jhoogenboom committed
515 516 517 518 519 520
                values = PAT_SPLIT.split(value)
                for value in values:
                    if PAT_SEQ_RAW.match(value) is None:
                        raise ValueError(
                            "Suffix sequence '%s' of marker %s is invalid" %
                            (value, marker))
521
                library["suffix"][marker] = values
jhoogenboom's avatar
jhoogenboom committed
522
                markers.add(marker)
523 524 525 526 527 528 529
            elif section_low == "genome_position":
                values = PAT_SPLIT.split(value)
                chromosome = PAT_CHROMOSOME.match(values[0])
                if chromosome is None:
                    raise ValueError(
                        "Invalid chromosome '%s' for marker %s." %
                        (values[0], marker))
530 531 532 533 534 535 536 537 538 539 540 541 542
                pos = [chromosome.group(1)]
                for i in range(1, len(values)):
                    try:
                        pos.append(int(values[i]))
                    except:
                        raise ValueError(
                            "Position '%s' of marker %s is not a valid integer"
                            % (values[i], marker))
                    if not i % 2 and pos[-2] >= pos[-1]:
                        raise ValueError(
                            "End position %i of marker %s must be higher than "
                            "corresponding start position %i" %
                            (pos[-1], marker, pos[-2]))
Hoogenboom, Jerry's avatar
Hoogenboom, Jerry committed
543 544
                if len(values) == 1:
                    pos.append(1)
545
                library["genome_position"][marker] = tuple(pos)
546
                markers.add(marker)
547
            elif section_low == "length_adjust":
jhoogenboom's avatar
jhoogenboom committed
548 549 550 551 552 553 554 555
                try:
                    value = int(value)
                except:
                    raise ValueError(
                        "Length adjustment '%s' of marker %s is not a valid "
                        "integer" % (value, marker))
                library["length_adjust"][marker] = value
                markers.add(marker)
556
            elif section_low == "block_length":
jhoogenboom's avatar
jhoogenboom committed
557 558 559 560 561 562 563 564
                try:
                    value = int(value)
                except:
                    raise ValueError(
                        "Block length '%s' of marker %s is not a valid integer"
                        % (value, marker))
                library["block_length"][marker] = value
                markers.add(marker)
565 566 567 568 569 570 571 572 573 574
            elif section_low == "max_expected_copies":
                try:
                    value = int(value)
                except:
                    raise ValueError(
                        "Maximum number of expected copies '%s' of marker %s "
                        "is not a valid integer" % (value, marker))
                library["max_expected_copies"][marker] = value
                markers.add(marker)
            elif section_low == "aliases":
jhoogenboom's avatar
jhoogenboom committed
575 576 577 578 579
                values = PAT_SPLIT.split(value)
                if len(values) != 3:
                    raise ValueError("Alias %s does not have 3 values, but %i"
                                     % (marker, len(values)))
                if PAT_SEQ_RAW.match(values[1]) is None:
580 581 582 583 584 585 586
                    raise ValueError(
                        "Alias sequence '%s' of alias %s is invalid" %
                        (values[1], marker))
                if PAT_ALIAS.match(values[2]) is None:
                    raise ValueError(
                        "Allele name '%s' of alias %s is invalid" %
                        (values[2], marker))
jhoogenboom's avatar
jhoogenboom committed
587 588 589 590 591 592
                library["aliases"][marker] = {
                    "marker": values[0],
                    "sequence": values[1],
                    "name": values[2]
                }
                markers.add(marker)
593
            elif section_low == "repeat":
594 595 596 597
                if marker in library["nostr_reference"]:
                    raise ValueError(
                        "Marker %s was encountered in both [repeat] and "
                        "[no_repeat] sections" % marker)
jhoogenboom's avatar
jhoogenboom committed
598 599 600 601 602 603
                if PAT_STR_DEF.match(value) is None:
                    raise ValueError(
                        "STR definition '%s' of marker %s is invalid" %
                        (value, marker))
                library["regex"][marker] = value
                markers.add(marker)
604 605 606 607 608 609 610 611 612 613 614 615 616 617 618
            elif section_low == "no_repeat":
                if marker in library["regex"]:
                    raise ValueError(
                        "Marker %s was encountered in both [repeat] and "
                        "[no_repeat] sections" % marker)
                if marker in library["prefix"] or marker in library["suffix"]:
                    raise ValueError(
                        "A prefix or suffix was defined for non-STR marker %s"
                        % marker)
                if PAT_SEQ_RAW.match(value) is None:
                    raise ValueError(
                        "Reference sequence '%s' of marker %s is invalid" %
                        (value, marker))
                library["nostr_reference"][marker] = value
                markers.add(marker)
619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637
            elif section_low == "expected_allele_length":
                values = PAT_SPLIT.split(value)
                try:
                    min_length = int(values[0])
                except:
                    raise ValueError(
                        "Minimum expected allele length '%s' of marker %s "
                        "is not a valid integer" % (values[0], marker))
                if len(values) > 1:
                    try:
                        max_length = int(values[1])
                    except:
                        raise ValueError(
                            "Maximum expected allele length '%s' of marker %s "
                            "is not a valid integer" % (values[1], marker))
                else:
                    max_length = sys.maxint
                library["expected_length"][marker] = (min_length, max_length)
                markers.add(marker)
638 639 640 641 642 643 644 645

    # Sanity check: prohibit prefix/suffix for aliases of non-STRs.
    for alias in library["aliases"]:
        if library["aliases"][alias]["marker"] in library["nostr_reference"] \
                and (alias in library["prefix"] or alias in library["suffix"]):
            raise ValueError(
                "A prefix or suffix was defined for alias %s of non-STR "
                "marker %s" % (alias, library["aliases"][alias]["marker"]))
jhoogenboom's avatar
jhoogenboom committed
646

647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
    # Sanity check: end position of marker should reflect ref length.
    for marker in library["genome_position"]:
        if marker not in library["nostr_reference"]:
            continue
        pos = library["genome_position"][marker]
        reflength = len(library["nostr_reference"][marker])
        length = 0
        for i in range(2, len(pos), 2):
            length += pos[i] - pos[i-1] + 1
        if reflength < length or (len(pos) % 2 and reflength != length):
            raise ValueError(
                "Length of reference sequence of marker %s is %i bases, but "
                "genome positions add up to %i bases" %
                (marker, reflength, length))

jhoogenboom's avatar
jhoogenboom committed
662 663 664 665
    # Compile regular expressions.
    # NOTE: The libconvert tool expects "(seq){num,num}" blocks ONLY!
    for marker in markers:
        parts = []
666
        blocksm = []
jhoogenboom's avatar
jhoogenboom committed
667
        if marker in library["prefix"]:
jhoogenboom's avatar
jhoogenboom committed
668
            parts += ("(%s){0,1}" % x for x in library["prefix"][marker])
jhoogenboom's avatar
jhoogenboom committed
669
        if marker in library["aliases"]:
670
            blocksm.append((library["aliases"][marker]["sequence"], 0, 1))
jhoogenboom's avatar
jhoogenboom committed
671
        elif marker in library["regex"]:
672 673 674
            blocksm += [(block[0], int(block[1]), int(block[2])) for block in
                        PAT_STR_DEF_BLOCK.findall(library["regex"][marker])]
        parts += ["(%s){%s,%s}" % x for x in blocksm]
jhoogenboom's avatar
jhoogenboom committed
675
        if marker in library["suffix"]:
jhoogenboom's avatar
jhoogenboom committed
676
            parts += ("(%s){0,1}" % x for x in library["suffix"][marker])
jhoogenboom's avatar
jhoogenboom committed
677 678 679
        if parts:
            library["regex"][marker] = re.compile(
                "".join(["^"] + parts + ["$"]))
680 681
        if blocksm:
            library["blocks_middle"][marker] = blocksm
jhoogenboom's avatar
jhoogenboom committed
682 683
    return library
#parse_library_ini
jhoogenboom's avatar
jhoogenboom committed
684 685


686
def load_profiles(profilefile, library=None):
jhoogenboom's avatar
jhoogenboom committed
687
    column_names = profilefile.readline().rstrip("\r\n").split("\t")
688 689 690
    (colid_marker, colid_allele, colid_sequence, colid_fmean, colid_rmean,
     colid_tool) = get_column_ids(column_names, "marker", "allele", "sequence",
        "fmean", "rmean", "tool")
jhoogenboom's avatar
jhoogenboom committed
691 692 693 694 695 696 697 698 699 700 701 702 703

    profiles = {}
    for line in profilefile:
        line = line.rstrip("\r\n").split("\t")
        if line == [""]:
            continue
        marker = line[colid_marker]
        if marker not in profiles:
            profiles[marker] = {
                "m": set(),  # To be replaced by its length below.
                "n": set(),  # To be replaced by its length below.
                "seqs": [],
                "forward": {},  # To be replaced by a list below.
704 705
                "reverse": {},  # To be replaced by a list below.
                "tool": {}  # To be replaced by a list below.
jhoogenboom's avatar
jhoogenboom committed
706 707 708 709 710 711 712 713 714 715 716 717
                }
        allele = ensure_sequence_format(line[colid_allele], "raw",
            library=library, marker=marker)
        sequence = ensure_sequence_format(line[colid_sequence], "raw",
            library=library, marker=marker)
        if (allele, sequence) in profiles[marker]["forward"]:
            raise ValueError(
                "Invalid background noise profiles file: encountered "
                "multiple values for marker '%s' allele '%s' sequence '%s'" %
                (marker, allele, sequence))
        profiles[marker]["forward"][allele,sequence] = float(line[colid_fmean])
        profiles[marker]["reverse"][allele,sequence] = float(line[colid_rmean])
718
        profiles[marker]["tool"][allele, sequence] = line[colid_tool]
jhoogenboom's avatar
jhoogenboom committed
719 720 721 722 723 724 725 726 727 728
        profiles[marker]["m"].update((allele, sequence))
        profiles[marker]["n"].add(allele)

    # Check completeness and reorder true alleles.
    for marker in profiles:
        profiles[marker]["seqs"] = list(profiles[marker]["n"]) + \
            list(profiles[marker]["m"]-profiles[marker]["n"])
        profiles[marker]["n"] = len(profiles[marker]["n"])
        profiles[marker]["m"] = len(profiles[marker]["m"])
        newprofiles = {"forward": [], "reverse": []}
729
        tools = []
jhoogenboom's avatar
jhoogenboom committed
730 731 732 733
        for i in range(profiles[marker]["n"]):
            allele = profiles[marker]["seqs"][i]
            for direction in newprofiles:
                newprofiles[direction].append([0] * profiles[marker]["m"])
734
            tools.append([""] * profiles[marker]["m"])
jhoogenboom's avatar
jhoogenboom committed
735 736 737 738 739 740
            for j in range(profiles[marker]["m"]):
                sequence = profiles[marker]["seqs"][j]
                if (allele, sequence) in profiles[marker]["forward"]:
                    for direction in newprofiles:
                        newprofiles[direction][i][j] = \
                            profiles[marker][direction][allele, sequence]
741
                    tools[i][j] = profiles[marker]["tool"][allele, sequence]
jhoogenboom's avatar
jhoogenboom committed
742 743
        profiles[marker]["forward"] = newprofiles["forward"]
        profiles[marker]["reverse"] = newprofiles["reverse"]
744
        profiles[marker]["tool"] = tools
jhoogenboom's avatar
jhoogenboom committed
745 746

    return profiles
747
#load_profiles
748 749


750
def pattern_longest_match(pattern, subject):
751
    """Return the longest match of the pattern in the subject string."""
752 753 754 755 756
    # FIXME, this function tries only one match at each position in the
    # sequence, which is not neccessarily the longest match at that
    # position. For now, we'll search the reverse sequence as well.
    # Re-implement this without regular expressions to test all options.
    reverse = False
757 758
    match = None
    pos = 0
759
    pat = re.compile("".join("(%s){%i,%i}" % x for x in pattern))
760
    while pos < len(subject):
761
        m = pat.search(subject, pos)
762 763 764 765 766
        if m is None:
            break
        if match is None or m.end()-m.start() > match.end()-match.start():
            match = m
        pos = m.start() + 1
767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815

    # Try to find a longer match from the other end.
    if match is not None:
        subject = reverse_complement(subject)
        pos = 0
        pat = re.compile("".join(
            "(%s){%i,%i}" % (reverse_complement(x[0]), x[1], x[2])
            for x in reversed(pattern)))
        while pos < len(subject):
            m = pat.search(subject, pos)
            if m is None:
                break
            if m.end()-m.start() > match.end()-match.start():
                match = m
                reverse = True
            pos = m.start() + 1

    # Extract the blocks from the match.
    match = [] if match is None or not match.group() else reduce(
        lambda x, i: (
            x[0] + [match.group(i)]*((match.end(i)-x[1])/len(match.group(i))),
            match.end(i)) if match.group(i) else x,
        range(1, match.lastindex+1), ([], match.start()))[0]

    # Return the match in the same sequence orientation as the input.
    return map(reverse_complement, reversed(match)) if reverse else match
#pattern_longest_match


def pattern_longest_match_veryslow(pattern, subject):
    """Return the longest match of the pattern in the subject string."""
    longest = 0
    the_match = []
    # Generate all possible matching sequences for this pattern.
    #print("Finding match of pattern %r to sequence %s" % (pattern, subject))
    for matching_blocks in itertools.product(*(
            [[block[0]]*i for i in range(block[1], block[2]+1)]
            for block in pattern)):
        matching = itertools.chain.from_iterable(matching_blocks)
        matching_seq = "".join(matching)
        matching_len = len(matching_seq)
        if matching_len <= longest:
            continue
        if matching_seq in subject:
            longest = matching_len
            the_match = matching
    #print("Found match covering %i/%i bases" % (longest, len(subject)))
    return the_match
#pattern_longest_match_veryslow
816 817


jhoogenboom's avatar
jhoogenboom committed
818 819 820
def detect_sequence_format(seq):
    """Return format of seq.  One of 'raw', 'tssv', or 'allelename'."""
    if not seq:
jhoogenboom's avatar
jhoogenboom committed
821
        raise ValueError("Empty sequence")
822
    if seq in SEQ_SPECIAL_VALUES:
823 824
        # Special case.
        return False
jhoogenboom's avatar
jhoogenboom committed
825 826 827 828
    if PAT_SEQ_RAW.match(seq):
        return 'raw'
    if PAT_SEQ_TSSV.match(seq):
        return 'tssv'
829 830
    if PAT_SEQ_ALLELENAME_STR.match(seq) or PAT_SEQ_ALLELENAME_MT.match(seq) \
            or PAT_SEQ_ALLELENAME_SNP.match(seq):
jhoogenboom's avatar
jhoogenboom committed
831
        return 'allelename'
jhoogenboom's avatar
jhoogenboom committed
832
    raise ValueError("Unrecognised sequence format")
jhoogenboom's avatar
jhoogenboom committed
833 834 835
#detect_sequence_format


jhoogenboom's avatar
jhoogenboom committed
836
def convert_sequence_tssv_raw(seq):
jhoogenboom's avatar
jhoogenboom committed
837 838
    return "".join(block[0] * int(block[1])
                   for block in PAT_TSSV_BLOCK.findall(seq))
jhoogenboom's avatar
jhoogenboom committed
839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854
#convert_sequence_tssv_raw


def convert_sequence_raw_tssv(seq, library, marker, return_alias=False):
    # Try to match this marker's pattern, or any of its aliases.
    match = None
    if "aliases" in library:
        for alias in library["aliases"]:
            if (library["aliases"][alias]["marker"] == marker and
                    alias in library["regex"]):
                match = library["regex"][alias].match(seq)
                if match is not None:
                    marker = alias
                    break
    if match is None and marker in library["regex"]:
        match = library["regex"][marker].match(seq)
855
    if match is not None:
jhoogenboom's avatar
jhoogenboom committed
856 857
        parts = ((match.group(i), match.end(i)) for i in range(1, 1 if
            match.lastindex is None else match.lastindex+1) if match.group(i))
jhoogenboom's avatar
jhoogenboom committed
858

859
    # Use heuristics if the sequence does not match the pattern.
jhoogenboom's avatar
jhoogenboom committed
860
    else:
861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877

        # Find explictily provided prefix and/or suffix if present.
        pre_suf = ["", ""]
        if "prefix" in library and marker in library["prefix"]:
            for prefix in library["prefix"][marker]:
                if seq.startswith(prefix):
                    pre_suf[0] = prefix
                    seq = seq[len(prefix):]
                    break
        if "suffix" in library and marker in library["suffix"]:
            for suffix in library["suffix"][marker]:
                if seq.endswith(suffix):
                    pre_suf[1] = suffix
                    seq = seq[:-len(suffix)]
                    break

        # Find longest match of middle pattern.
878
        middle = [(seq, len(pre_suf[0])+len(seq))] if seq else []
879 880 881 882
        if middle and marker in library["blocks_middle"]:
            match = pattern_longest_match(library["blocks_middle"][marker],seq)
            matched = "".join(match)
            if matched:
883 884 885 886 887 888

                # If this allele does not match the prefix of this
                # marker, but the canonical prefix of the marker ends
                # with the same sequence as the start of our match, we
                # move that portion of the match into the prefix.
                # Then, we do the same thing with the suffix.
889 890 891 892 893
                middle = match
                match_start = seq.index(matched)
                match_end = match_start + len(matched)
                start = match_start
                end = match_end
894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917
                modified = False
                if (not pre_suf[0] and "prefix" in library
                        and marker in library["prefix"]):
                    ref = library["prefix"][marker][0]
                    i = min(len(ref), len(matched))
                    while i > 0:
                        if ref.endswith(matched[:i]):
                            start += i
                            matched = matched[i:]
                            modified = True
                            break
                        i -= 1
                if (not pre_suf[1] and "suffix" in library
                        and marker in library["suffix"]):
                    ref = library["suffix"][marker][0]
                    i = min(len(ref), len(matched))
                    while i > 0:
                        if ref.startswith(matched[-i:]):
                            end -= i
                            matched = matched[:-i]
                            modified = True
                            break
                        i -= 1
                if modified:
918 919
                    from_start = start - match_start
                    from_end = match_end - end
920 921 922 923 924 925 926 927 928 929 930 931 932 933
                    while from_start:
                        if from_start < len(middle[0]):
                            middle[0] = middle[0][from_start:]
                            break
                        else:
                            from_start -= len(middle[0])
                            middle = middle[1:]
                    while from_end:
                        if from_end < len(middle[-1]):
                            middle[-1] = middle[-1][:-from_end]
                            break
                        else:
                            from_end -= len(middle[-1])
                            middle = middle[:-1]
934 935 936 937 938 939
                if middle:
                    middle = reduce(
                        lambda x, y: (x[:-1] if x[-1][0] == y else x) +
                            [(y, x[-1][1]+len(y))], middle[1:],
                            [(middle[0],
                              start+len(middle[0])+len(pre_suf[0]))])
940 941 942 943

                pre_suf[0] += seq[:start]
                pre_suf[1] = seq[end:] + pre_suf[1]
                seq = matched
944 945 946 947 948 949 950 951 952 953 954 955 956 957 958

        # Now construct parts.
        parts = []
        if pre_suf[0]:
            parts.append((pre_suf[0], len(pre_suf[0])))
        parts += middle
        if pre_suf[1]:
            parts.append((pre_suf[1], sum(map(len,pre_suf))+len(seq)))

    seq = reduce(
        lambda a, b: (a[0] + "%s(%i)" % (b[0], (b[1]-a[1])/len(b[0])), b[1]),
        reduce(
            lambda x, y: x[:-1] + [y] if x[-1][0] == y[0] else x + [y],
            parts,
            [("", 0)]))[0]
jhoogenboom's avatar
jhoogenboom committed
959 960 961 962 963 964
    return (seq, marker) if return_alias else seq
#convert_sequence_raw_tssv


def convert_sequence_allelename_tssv(seq, library, marker):
    # Check whether there is an alias for this sequence.
965
    alias_of = None
jhoogenboom's avatar
jhoogenboom committed
966 967 968 969 970
    if "aliases" in library:
        for alias in library["aliases"]:
            if library["aliases"][alias]["marker"] == marker and (
                    seq == library["aliases"][alias]["name"] or
                    seq.startswith(library["aliases"][alias]["name"] + "_")):
971
                alias_of = marker
jhoogenboom's avatar
jhoogenboom committed
972 973 974 975 976 977 978
                marker = alias
                seq = "".join([
                    "0_",
                    library["aliases"][alias]["sequence"] + "[1]",
                    seq[len(library["aliases"][alias]["name"]):]])
                break

979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001
    nameformat = None
    if PAT_SEQ_ALLELENAME_MT.match(seq) is not None:
        nameformat = "MtDNA"
    elif PAT_SEQ_ALLELENAME_SNP.match(seq) is not None:
        nameformat = "SNP"
    if nameformat is not None:
        # MtDNA and SNP markers.
        try:
            reference = library["nostr_reference"][marker]
        except KeyError:
            raise ValueError(
                "%s allele '%s' found for marker %s, but "
                "no reference sequence was found in the library" %
                (nameformat, seq, marker))
        if seq == "REF":
            return reference + "(1)"
        return mutate_sequence(reference, seq.split(),
            library["genome_position"].get(marker,
                ("M" if nameformat == "MtDNA" else "", 1))) + "(1)"

    # Note: aliases of mtDNA and SNP markers end up here as well.
    # It should NOT look like an alias now, however.
    match = PAT_SEQ_ALLELENAME_STR.match(seq)
1002 1003 1004
    if match is None or match.group(1) is not None:
        raise ValueError("Invalid allele name '%s' encountered!" % seq)

jhoogenboom's avatar
jhoogenboom committed
1005 1006 1007
    allele = seq.split("_")

    # Get and mutate prefix and suffix.
1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019
    prefix = ""
    suffix = ""
    if "prefix" in library:
        if marker in library["prefix"]:
            prefix = library["prefix"][marker][0]
        elif alias_of is not None and alias_of in library["prefix"]:
            prefix = library["prefix"][alias_of][0]
    if "suffix" in library:
        if marker in library["suffix"]:
            suffix = library["suffix"][marker][0]
        elif alias_of is not None and alias_of in library["suffix"]:
            suffix = library["suffix"][alias_of][0]
jhoogenboom's avatar
jhoogenboom committed
1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045
    variants = [[], []]
    for variant in allele[2:]:
        if variant[0] == "-":
            if not prefix:
                raise ValueError("Encountered prefix variant '%s', but marker "
                                 "'%s' has no prefix!" % (variant, marker))
            variants[0].append(variant)
        elif variant[0] == "+":
            if not suffix:
                raise ValueError("Encountered suffix variant '%s', but marker "
                                 "'%s' has no suffix!" % (variant, marker))
            variants[1].append(variant)
        else:
            raise ValueError("Unrecognised variant '%s'" % variant)
    if variants[0]:
        prefix = mutate_sequence(prefix, variants[0])
    if variants[1]:
        suffix = mutate_sequence(suffix, variants[1])

    blocks = []
    if prefix:
        blocks.append((prefix, 1))
    for block in PAT_ALLELENAME_BLOCK.findall(allele[1]):
        blocks.append((block[0], int(block[1])))
    if suffix:
        blocks.append((suffix, 1))
jhoogenboom's avatar
jhoogenboom committed
1046
    return "".join("%s(%i)" % block for block in blocks)
jhoogenboom's avatar
jhoogenboom committed
1047 1048 1049 1050 1051 1052
#convert_sequence_allelename_tssv


def convert_sequence_raw_allelename(seq, library, marker):
    # We actually convert raw->allelename via TSSV format.
    seq, alias = convert_sequence_raw_tssv(seq, library, marker, True)
1053 1054
    blocks = PAT_TSSV_BLOCK.findall(seq)

1055
    if "nostr_reference" in library and marker in library["nostr_reference"]:
1056 1057 1058
        # Handle non-STR markers here.
        if alias != marker:
            return library["aliases"][alias]["name"]
1059 1060 1061
        if not blocks:
            # Oh dear, empty sequence... Primer dimer?
            blocks = (("",),)
1062 1063 1064 1065 1066 1067
        if library["nostr_reference"][marker] == blocks[0][0]:
            return "REF"
        return " ".join(
            call_variants(library["nostr_reference"][marker], blocks[0][0],
                library["genome_position"].get(marker, "suffix")))

1068 1069 1070
    # Find prefix and suffix.
    prefix = suffix = this_prefix = this_suffix = ""
    remaining_blocks = len(blocks)
1071 1072 1073 1074 1075
    if "prefix" in library:
        if alias in library["prefix"]:
            prefix = library["prefix"][alias][0]
        elif marker in library["prefix"]:
            prefix = library["prefix"][marker][0]
jhoogenboom's avatar
jhoogenboom committed
1076
        if prefix and remaining_blocks > 0 and blocks[0][1] == "1":
1077
            remaining_blocks -= 1
1078 1079 1080 1081 1082
    if "suffix" in library:
        if alias in library["suffix"]:
            suffix = library["suffix"][alias][0]
        elif marker in library["suffix"]:
            suffix = library["suffix"][marker][0]
jhoogenboom's avatar
jhoogenboom committed
1083
        if suffix and remaining_blocks > 0 and blocks[-1][1] == "1":
1084
            remaining_blocks -= 1
1085 1086 1087 1088 1089 1090
    if remaining_blocks > 0 and prefix and blocks[0][1] == "1":
        this_prefix = blocks[0][0]
        blocks = blocks[1:]
    if remaining_blocks > 0 and suffix and blocks[-1][1] == "1":
        this_suffix = blocks[-1][0]
        blocks = blocks[:-1]
jhoogenboom's avatar
jhoogenboom committed
1091 1092 1093 1094

    # Generate prefix/suffix variants.
    length = 0
    variants = []
1095
    if prefix != this_prefix:
1096
        variants += call_variants(prefix, this_prefix, "prefix")
1097 1098
        length += len(this_prefix) - len(prefix)
    if suffix != this_suffix:
1099
        variants += call_variants(suffix, this_suffix, "suffix")
1100