lib.py 66.1 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)
175 176 177 178 179
    elif location[0] == "M":
        MATCH_SCORE = 1
        MISMATCH_SCORE = -1
        GAP_OPEN_SCORE = -2
        GAP_EXTEND_SCORE = -1
jhoogenboom's avatar
jhoogenboom committed
180 181 182 183 184 185 186 187 188 189 190

    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]
191
                matrix_direction[i] = A_HORZ_E | (A_HORZ_O if x == 1 else 0)
jhoogenboom's avatar
jhoogenboom committed
192 193 194 195
            elif y != 0:
                # Left column.
                matrix_gap2[i] = GAP_OPEN_SCORE + GAP_EXTEND_SCORE * (y - 1)
                matrix_match[i] = matrix_gap2[i]
196 197 198 199
                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
200 201 202 203 204 205 206
            continue

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

207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
        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
222 223 224
            matrix_match[i-1-row_offset] + match,
            matrix_gap1[i],
            matrix_gap2[i])
225 226 227 228 229 230 231
        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
232 233 234 235 236 237 238 239 240 241 242

    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]))
243 244 245 246 247 248 249 250 251 252
        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
253 254 255 256 257 258 259


    # Backtracking.
    variants = []
    variant_template = 0
    variant_sequence = 0
    i = len(matrix_match) - 1
260
    in_matrix = M_MATCH  # May change before first step.
jhoogenboom's avatar
jhoogenboom committed
261 262 263
    while i >= 0:
        x = i % row_offset
        y = i / row_offset
264 265 266 267 268 269 270 271 272 273 274 275 276 277
        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
278 279
            # Go horizontally.  Deletion.
            variant_template += 1
280 281 282
            if matrix_direction[i] & A_HORZ_O:
                # End of gap, go diagonally after this.
                in_matrix = M_MATCH
jhoogenboom's avatar
jhoogenboom committed
283 284 285
            i -= 1
            continue

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

        # Go diagonally.  Either match or mismatch.
296 297 298 299 300 301
        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
302 303
            # Match.  Flush variants.
            if variant_template or variant_sequence:
304 305 306 307 308 309
                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 "",#"-",
310 311
                            get_genome_pos(
                                location, x + min(j, variant_template-1)),
312 313
                            ".%i" % (j-variant_template+1)
                                if j >= variant_template else "",
314
                            sequence[y+j] if j < variant_sequence else "del"))
315
                elif variant_template == 0:
jhoogenboom's avatar
jhoogenboom committed
316
                    # Insertions: "-131.1->C" instead of "-130->C".
317
                    variants.append(variant_format % (
318
                        get_genome_pos(location, x - 1),
319
                        ".1-",
jhoogenboom's avatar
jhoogenboom committed
320 321
                        sequence[y:y+variant_sequence]))
                else:
322
                    variants.append(variant_format % (
323
                        get_genome_pos(location, x),
jhoogenboom's avatar
jhoogenboom committed
324 325 326 327 328 329
                        template[x:x+variant_template],
                        sequence[y:y+variant_sequence] or "-"))
                variant_template = 0
                variant_sequence = 0
        i -= 1 + row_offset

330 331
    # Variants were called from right to left.  Reverse their order.
    if location[0] != "prefix":
jhoogenboom's avatar
jhoogenboom committed
332 333 334 335
        variants.reverse()

    # Store the result in the cache.
    if cache:
336
        call_variants.cache[template, sequence, cache_key] = variants
jhoogenboom's avatar
jhoogenboom committed
337 338 339 340 341
    return variants
#call_variants
call_variants.cache = {}


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

    seq = [[]] + [[base] for base in seq]
jhoogenboom's avatar
jhoogenboom committed
355
    for variant in variants:
356
        vm = pattern.match(variant)
jhoogenboom's avatar
jhoogenboom committed
357 358
        if vm is None:
            raise ValueError("Unrecognised variant '%s'" % variant)
359 360 361 362
        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
363 364
        if old == "-":
            old = ""
365
        if new == "-" or new == "del":
jhoogenboom's avatar
jhoogenboom committed
366 367
            new = ""
        if pos < 0:
368
            pos += len(seq)
369 370
        pos = get_genome_pos(location, pos, True)
        if pos < 0 or (pos == 0 and not ins) or pos >= len(seq):
371
            raise ValueError(
372 373
                "Position of variant '%s' is outside sequence range" %
                    (variant))
374 375 376 377 378 379 380 381 382 383 384 385 386 387 388
        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
389
        else:
390 391 392 393 394
            # 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
395 396 397
#mutate_sequence


398
def parse_library(libfile, stream=False):
jhoogenboom's avatar
jhoogenboom committed
399
    try:
400
        if not stream:
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415
            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:
416 417
            libfile.close()
        return library
418 419
    except ValueError as err:
        raise argparse.ArgumentTypeError(err)
jhoogenboom's avatar
jhoogenboom committed
420 421 422 423 424 425 426 427 428 429 430 431 432 433
#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": {},
434
      "regex": {},
435
      "blocks_middle": {}
jhoogenboom's avatar
jhoogenboom committed
436 437
    }
    for line in handle:
jhoogenboom's avatar
jhoogenboom committed
438
        line = [x.strip() for x in line.rstrip("\r\n").split("\t")]
jhoogenboom's avatar
jhoogenboom committed
439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455
        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]
456 457 458
        library["blocks_middle"][marker] = [
            (block[0], int(block[1]), int(block[2])) for block in
                PAT_STR_DEF_BLOCK.findall(line[3])]
459
        library["regex"][marker] = re.compile(
460 461 462
            "".join(("^", "".join(
                "(%s){%s,%s}" % x for x in PAT_STR_DEF_BLOCK.findall(line[3])),
                "$")))
jhoogenboom's avatar
jhoogenboom committed
463 464 465 466 467 468 469 470 471 472
    return library
#parse_library_tsv


def parse_library_ini(handle):
    library = {
      "flanks": {},
      "prefix": {},
      "suffix": {},
      "regex": {},
473
      "blocks_middle": {},
474 475
      "nostr_reference": {},
      "genome_position": {},
476 477
      "length_adjust": {},
      "block_length": {},
478
      "max_expected_copies": {},
479
      "expected_length": {},
jhoogenboom's avatar
jhoogenboom committed
480 481 482 483 484 485 486 487 488 489
      "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)
490 491
            section_low = section.lower()
            if section_low == "flanks":
jhoogenboom's avatar
jhoogenboom committed
492 493 494 495 496 497 498 499 500 501 502 503
                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)
504
            elif section_low == "prefix":
505 506 507
                if marker in library["nostr_reference"]:
                    raise ValueError(
                        "A prefix was defined for non-STR marker %s" % marker)
jhoogenboom's avatar
jhoogenboom committed
508 509 510 511 512 513
                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))
514
                library["prefix"][marker] = values
jhoogenboom's avatar
jhoogenboom committed
515
                markers.add(marker)
516
            elif section_low == "suffix":
517 518 519
                if marker in library["nostr_reference"]:
                    raise ValueError(
                        "A suffix was defined for non-STR marker %s" % marker)
jhoogenboom's avatar
jhoogenboom committed
520 521 522 523 524 525
                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))
526
                library["suffix"][marker] = values
jhoogenboom's avatar
jhoogenboom committed
527
                markers.add(marker)
528 529 530 531 532 533 534
            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))
535 536 537 538 539 540 541 542 543 544 545 546 547
                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
548 549
                if len(values) == 1:
                    pos.append(1)
550
                library["genome_position"][marker] = tuple(pos)
551
                markers.add(marker)
552
            elif section_low == "length_adjust":
jhoogenboom's avatar
jhoogenboom committed
553 554 555 556 557 558 559 560
                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)
561
            elif section_low == "block_length":
jhoogenboom's avatar
jhoogenboom committed
562 563 564 565 566 567 568 569
                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)
570 571 572 573 574 575 576 577 578 579
            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
580 581 582 583 584
                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:
585 586 587 588 589 590 591
                    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
592 593 594 595 596 597
                library["aliases"][marker] = {
                    "marker": values[0],
                    "sequence": values[1],
                    "name": values[2]
                }
                markers.add(marker)
598
            elif section_low == "repeat":
599 600 601 602
                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
603 604 605 606 607 608
                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)
609 610 611 612 613 614 615 616 617 618 619 620 621 622 623
            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)
624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642
            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)
643 644 645 646 647 648 649 650

    # 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
651

652 653 654 655 656 657 658 659 660 661 662 663 664 665 666
    # 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
667 668 669 670
    # Compile regular expressions.
    # NOTE: The libconvert tool expects "(seq){num,num}" blocks ONLY!
    for marker in markers:
        parts = []
671
        blocksm = []
jhoogenboom's avatar
jhoogenboom committed
672
        if marker in library["prefix"]:
jhoogenboom's avatar
jhoogenboom committed
673
            parts += ("(%s){0,1}" % x for x in library["prefix"][marker])
jhoogenboom's avatar
jhoogenboom committed
674
        if marker in library["aliases"]:
675
            blocksm.append((library["aliases"][marker]["sequence"], 0, 1))
jhoogenboom's avatar
jhoogenboom committed
676
        elif marker in library["regex"]:
677 678 679
            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
680
        if marker in library["suffix"]:
jhoogenboom's avatar
jhoogenboom committed
681
            parts += ("(%s){0,1}" % x for x in library["suffix"][marker])
jhoogenboom's avatar
jhoogenboom committed
682 683 684
        if parts:
            library["regex"][marker] = re.compile(
                "".join(["^"] + parts + ["$"]))
685 686
        if blocksm:
            library["blocks_middle"][marker] = blocksm
jhoogenboom's avatar
jhoogenboom committed
687 688
    return library
#parse_library_ini
jhoogenboom's avatar
jhoogenboom committed
689 690


691
def load_profiles(profilefile, library=None):
jhoogenboom's avatar
jhoogenboom committed
692
    column_names = profilefile.readline().rstrip("\r\n").split("\t")
693 694 695
    (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
696 697 698 699 700 701 702 703 704 705 706 707 708

    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.
709 710
                "reverse": {},  # To be replaced by a list below.
                "tool": {}  # To be replaced by a list below.
jhoogenboom's avatar
jhoogenboom committed
711 712 713 714 715 716 717 718 719 720 721 722
                }
        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])
723
        profiles[marker]["tool"][allele, sequence] = line[colid_tool]
jhoogenboom's avatar
jhoogenboom committed
724 725 726 727 728 729 730 731 732 733
        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": []}
734
        tools = []
jhoogenboom's avatar
jhoogenboom committed
735 736 737 738
        for i in range(profiles[marker]["n"]):
            allele = profiles[marker]["seqs"][i]
            for direction in newprofiles:
                newprofiles[direction].append([0] * profiles[marker]["m"])
739
            tools.append([""] * profiles[marker]["m"])
jhoogenboom's avatar
jhoogenboom committed
740 741 742 743 744 745
            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]
746
                    tools[i][j] = profiles[marker]["tool"][allele, sequence]
jhoogenboom's avatar
jhoogenboom committed
747 748
        profiles[marker]["forward"] = newprofiles["forward"]
        profiles[marker]["reverse"] = newprofiles["reverse"]
749
        profiles[marker]["tool"] = tools
jhoogenboom's avatar
jhoogenboom committed
750 751

    return profiles
752
#load_profiles
753 754


755
def pattern_longest_match(pattern, subject):
756
    """Return the longest match of the pattern in the subject string."""
757 758 759 760 761
    # 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
762 763
    match = None
    pos = 0
764
    pat = re.compile("".join("(%s){%i,%i}" % x for x in pattern))
765
    while pos < len(subject):
766
        m = pat.search(subject, pos)
767 768 769 770 771
        if m is None:
            break
        if match is None or m.end()-m.start() > match.end()-match.start():
            match = m
        pos = m.start() + 1
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 816 817 818 819 820

    # 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
821 822


jhoogenboom's avatar
jhoogenboom committed
823 824 825
def detect_sequence_format(seq):
    """Return format of seq.  One of 'raw', 'tssv', or 'allelename'."""
    if not seq:
jhoogenboom's avatar
jhoogenboom committed
826
        raise ValueError("Empty sequence")
827
    if seq in SEQ_SPECIAL_VALUES:
828 829
        # Special case.
        return False
jhoogenboom's avatar
jhoogenboom committed
830 831 832 833
    if PAT_SEQ_RAW.match(seq):
        return 'raw'
    if PAT_SEQ_TSSV.match(seq):
        return 'tssv'
834 835
    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
836
        return 'allelename'
jhoogenboom's avatar
jhoogenboom committed
837
    raise ValueError("Unrecognised sequence format")
jhoogenboom's avatar
jhoogenboom committed
838 839 840
#detect_sequence_format


jhoogenboom's avatar
jhoogenboom committed
841
def convert_sequence_tssv_raw(seq):
jhoogenboom's avatar
jhoogenboom committed
842 843
    return "".join(block[0] * int(block[1])
                   for block in PAT_TSSV_BLOCK.findall(seq))
jhoogenboom's avatar
jhoogenboom committed
844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859
#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)
860
    if match is not None:
jhoogenboom's avatar
jhoogenboom committed
861 862
        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
863

864
    # Use heuristics if the sequence does not match the pattern.
jhoogenboom's avatar
jhoogenboom committed
865
    else:
866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882

        # 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.
883
        middle = [(seq, len(pre_suf[0])+len(seq))] if seq else []
884 885 886 887
        if middle and marker in library["blocks_middle"]:
            match = pattern_longest_match(library["blocks_middle"][marker],seq)
            matched = "".join(match)
            if matched:
888 889 890 891 892 893

                # 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.
894 895 896 897 898
                middle = match
                match_start = seq.index(matched)
                match_end = match_start + len(matched)
                start = match_start
                end = match_end
899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922
                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:
923 924
                    from_start = start - match_start
                    from_end = match_end - end
925 926 927 928 929 930 931 932 933 934 935 936 937 938
                    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]
939 940 941 942 943 944
                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]))])
945 946 947 948

                pre_suf[0] += seq[:start]
                pre_suf[1] = seq[end:] + pre_suf[1]
                seq = matched
949 950 951 952 953 954 955 956 957 958 959 960 961 962 963

        # 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
964 965 966 967 968 969
    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.
970
    alias_of = None
jhoogenboom's avatar
jhoogenboom committed
971 972 973 974 975
    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"] + "_")):
976
                alias_of = marker
jhoogenboom's avatar
jhoogenboom committed
977 978 979 980 981 982 983
                marker = alias
                seq = "".join([
                    "0_",
                    library["aliases"][alias]["sequence"] + "[1]",
                    seq[len(library["aliases"][alias]["name"]):]])
                break

984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006
    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)
1007 1008 1009
    if match is None or match.group(1) is not None:
        raise ValueError("Invalid allele name '%s' encountered!" % seq)

jhoogenboom's avatar
jhoogenboom committed
1010 1011 1012
    allele = seq.split("_")

    # Get and mutate prefix and suffix.
1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024
    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
1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050
    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
1051
    return "".join("%s(%i)" % block for block in blocks)
jhoogenboom's avatar
jhoogenboom committed
1052 1053 1054 1055 1056 1057
#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)
1058 1059
    blocks = PAT_TSSV_BLOCK.findall(seq)

1060
    if "nostr_reference" in library and marker in library["nostr_reference"]:
1061 1062 1063
        # Handle non-STR markers here.
        if alias != marker:
            return library["aliases"][alias]["name"]
1064 1065 1066
        if not blocks:
            # Oh dear, empty sequence... Primer dimer?
            blocks = (("",),)
1067 1068 1069 1070 1071 1072
        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")))

1073 1074 1075
    # Find prefix and suffix.
    prefix = suffix = this_prefix = this_suffix = ""
    remaining_blocks = len(blocks)
1076 1077 1078 1079 1080
    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
1081
        if prefix and remaining_blocks > 0 and blocks[0][1] == "1":
1082
            remaining_blocks -= 1
1083 1084 1085 1086 1087
    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
1088
        if suffix and remaining_blocks > 0 and blocks[-1][1] == "1":
1089
            remaining_blocks -= 1
1090 1091 1092 1093 1094 1095
    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
1096 1097 1098 1099

    # Generate prefix/suffix variants.
    length = 0
    variants = []
1100
    if prefix != this_prefix:
1101
        variants += call_variants(prefix, this_prefix, "prefix")
1102 1103
        length += len(this_prefix) - len(prefix)
    if suffix != this_suffix: