duplex: utils/sim.py annotate

author	nick
date	Thu, 02 Feb 2017 18:44:31 -0500
parents
children

rev	line source
18 e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	1 #!/usr/bin/env python
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	2 from __future__ import division
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	3 from __future__ import print_function
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	4 import re
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	5 import os
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	6 import sys
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	7 import copy
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	8 import numpy
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	9 import bisect
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	10 import random
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	11 import string
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	12 import numbers
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	13 import tempfile
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	14 import argparse
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	15 import subprocess
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	16 import fastqreader
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	17
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	18 REVCOMP_TABLE = string.maketrans('acgtrymkbdhvACGTRYMKBDHV', 'tgcayrkmvhdbTGCAYRKMVHDB')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	19 WGSIM_ID_REGEX = r'^(.+)_(\d+)_(\d+)_\d+:\d+:\d+_\d+:\d+:\d+_([0-9a-f]+)/[12]$'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	20 ARG_DEFAULTS = {'read_len':100, 'frag_len':400, 'n_frags':1000, 'out_format':'fasta',
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	21 'seq_error':0.001, 'pcr_error':0.001, 'cycles':25, 'indel_rate':0.15,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	22 'ext_rate':0.3, 'seed':None, 'invariant':'TGACT', 'bar_len':12, 'fastq_qual':'I'}
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	23 USAGE = "%(prog)s [options]"
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	24 DESCRIPTION = """Simulate a duplex sequencing experiment."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	25
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	26 RAW_DISTRIBUTION = (
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	27 # 0 1 2 3 4 5 6 7 8 9
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	28 # Low singletons, but then constant drop-off. From pML113 (see 2015-09-28 report).
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	29 # 0, 100, 36, 31, 27, 22, 17, 12, 7, 4.3,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	30 #2.4, 1.2, 0.6, 0.3, 0.2, 0.15, 0.1, 0.07, 0.05, 0.03,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	31 # High singletons, but then a second peak around 10. From Christine plasmid (2015-10-06 report).
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	32 # 0, 100, 5.24, 3.67, 3.50, 3.67, 3.85, 4.02, 4.11, 4.20,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	33 # 4.17, 4.10, 4.00, 3.85, 3.69, 3.55, 3.38, 3.15, 2.92, 2.62,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	34 # 2.27, 2.01, 1.74, 1.56, 1.38, 1.20, 1.02, 0.85,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	35 # Same as above, but low singletons, 2's, and 3's (rely on errors to fill out those).
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	36 0, 1, 2, 3, 3.50, 3.67, 3.85, 4.02, 4.11, 4.20,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	37 4.17, 4.10, 4.00, 3.85, 3.69, 3.55, 3.38, 3.15, 2.92, 2.62,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	38 2.27, 2.01, 1.74, 1.56, 1.38, 1.20, 1.02, 0.85,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	39 )
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	40
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	41
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	42 def main(argv):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	43
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	44 parser = argparse.ArgumentParser(description=DESCRIPTION)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	45 parser.set_defaults(**ARG_DEFAULTS)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	46
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	47 parser.add_argument('ref', metavar='ref.fa', nargs='?',
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	48 help='Reference sequence. Omit if giving --frag-file.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	49 parser.add_argument('out1', type=argparse.FileType('w'),
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	50 help='Write final mate 1 reads to this file.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	51 parser.add_argument('out2', type=argparse.FileType('w'),
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	52 help='Write final mate 2 reads to this file.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	53 parser.add_argument('-o', '--out-format', choices=('fastq', 'fasta'))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	54 parser.add_argument('--stdout', action='store_true',
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	55 help='Print interleaved output reads to stdout.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	56 parser.add_argument('-m', '--mutations', type=argparse.FileType('w'),
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	57 help='Write a log of the PCR and sequencing errors introduced to this file. Will overwrite any '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	58 'existing file at this path.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	59 parser.add_argument('-b', '--barcodes', type=argparse.FileType('w'),
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	60 help='Write a log of which barcodes were ligated to which fragments. Will overwrite any '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	61 'existing file at this path.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	62 parser.add_argument('--frag-file',
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	63 help='The path of the FASTQ file of fragments. If --ref is given, these will be generated with '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	64 'wgsim and kept (normally a temporary file is used, then deleted). Note: the file will be '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	65 'overwritten! If --ref is not given, then this should be a file of already generated '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	66 'fragments, and they will be used instead of generating new ones.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	67 parser.add_argument('-Q', '--fastq-qual',
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	68 help='The quality score to assign to all bases in FASTQ output. Give a character or PHRED '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	69 'score (integer). A PHRED score will be converted using the Sanger offset (33). Default: '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	70 '"%(default)s"')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	71 parser.add_argument('-S', '--seed', type=int,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	72 help='Random number generator seed. By default, a random, 32-bit seed will be generated and '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	73 'logged to stdout.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	74 params = parser.add_argument_group('simulation parameters')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	75 params.add_argument('-n', '--n-frags', type=int,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	76 help='The number of original fragment molecules to simulate. The final number of reads will be '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	77 'this multiplied by the average number of reads per family. If you provide fragments with '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	78 '--frag-file, the script will still only read in the number specified here. Default: '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	79 '%(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	80 params.add_argument('-r', '--read-len', type=int,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	81 help='Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	82 params.add_argument('-f', '--frag-len', type=int,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	83 help='Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	84 params.add_argument('-s', '--seq-error', type=float,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	85 help='Sequencing error rate per base (0-1 proportion, not percent). Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	86 params.add_argument('-p', '--pcr-error', type=float,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	87 help='PCR error rate per base (0-1 proportion, not percent). Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	88 params.add_argument('-c', '--cycles', type=int,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	89 help='Number of PCR cycles to simulate. Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	90 params.add_argument('-i', '--indel-rate', type=float,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	91 help='Fraction of errors which are indels. Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	92 params.add_argument('-E', '--extension-rate', dest='ext_rate', type=float,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	93 help='Probability an indel is extended. Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	94 params.add_argument('-B', '--bar-len', type=int,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	95 help='Length of the barcodes to generate. Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	96 params.add_argument('-I', '--invariant',
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	97 help='The invariant linker sequence between the barcode and sample sequence in each read. '
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	98 'Default: %(default)s')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	99
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	100 # Parse and interpret arguments.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	101 args = parser.parse_args(argv[1:])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	102 assert args.ref or args.frag_file, 'You must provide either a reference or fragments file.'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	103 if args.seed is None:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	104 seed = random.randint(0, 2**31-1)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	105 sys.stderr.write('seed: {}\n'.format(seed))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	106 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	107 seed = args.seed
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	108 random.seed(seed)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	109 if args.stdout:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	110 out1 = sys.stdout
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	111 out2 = sys.stdout
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	112 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	113 out1 = args.out1
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	114 out2 = args.out2
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	115 if isinstance(args.fastq_qual, numbers.Integral):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	116 assert args.fastq_qual >= 0, '--fastq-qual cannot be negative.'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	117 fastq_qual = chr(args.fastq_qual + 33)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	118 elif isinstance(args.fastq_qual, basestring):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	119 assert len(args.fastq_qual) == 1, '--fastq-qual cannot be more than a single character.'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	120 fastq_qual = args.fastq_qual
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	121 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	122 raise AssertionError('--fastq-qual must be a positive integer or single character.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	123 qual_line = fastq_qual * args.read_len
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	124
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	125 invariant_rc = get_revcomp(args.invariant)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	126
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	127 # Create a temporary director to do our work in. Then work inside a try so we can finally remove
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	128 # the directory no matter what exceptions are encountered.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	129 tmpfile = tempfile.NamedTemporaryFile(prefix='wgdsim.frags.')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	130 tmpfile.close()
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	131 try:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	132 # Step 1: Use wgsim to create fragments from the reference.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	133 if args.frag_file:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	134 frag_file = args.frag_file
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	135 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	136 frag_file = tmpfile.name
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	137 if args.ref and os.path.isfile(args.ref) and os.path.getsize(args.ref):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	138 #TODO: Check exit status
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	139 #TODO: Check for wgsim on the PATH.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	140 # Set error and mutation rates to 0 to just slice sequences out of the reference without
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	141 # modification.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	142 run_command('wgsim', '-e', '0', '-r', '0', '-d', '0', '-R', args.indel_rate, '-S', seed,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	143 '-N', args.n_frags, '-X', args.ext_rate, '-1', args.frag_len,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	144 args.ref, frag_file, os.devnull)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	145
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	146 # NOTE: Coordinates here are 0-based (0 is the first base in the sequence).
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	147 extended_dist = extend_dist(RAW_DISTRIBUTION)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	148 proportional_dist = compile_dist(extended_dist)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	149 n_frags = 0
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	150 for raw_fragment in fastqreader.FastqReadGenerator(frag_file):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	151 n_frags += 1
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	152 if n_frags > args.n_frags:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	153 break
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	154 chrom, id_num, start, stop = parse_read_id(raw_fragment.id)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	155 barcode1 = get_rand_seq(args.bar_len)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	156 barcode2 = get_rand_seq(args.bar_len)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	157 barcode2_rc = get_revcomp(barcode2)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	158 raw_frag_full = barcode1 + args.invariant + raw_fragment.seq + invariant_rc + barcode2
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	159
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	160 # Step 2: Determine how many reads to produce from each fragment.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	161 # - Use random.random() and divide the range 0-1 into segments of sizes proportional to
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	162 # the likelihood of each family size.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	163 # bisect.bisect() finds where an element belongs in a sorted list, returning the index.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	164 # proportional_dist is just such a sorted list, with values from 0 to 1.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	165 n_reads = bisect.bisect(proportional_dist, random.random())
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	166
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	167 # Step 3: Introduce PCR errors.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	168 # - Determine the mutations and their frequencies.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	169 # - Could get frequency from the cycle of PCR it occurs in.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	170 # - Important to have PCR errors shared between reads.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	171 # - For each read, determine which mutations it contains.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	172 # - Use random.random() < mut_freq.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	173 tree = get_good_pcr_tree(n_reads, args.cycles, 1000, max_diff=1)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	174 # Add errors to all children of original fragment.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	175 subtree1 = tree.get('child1')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	176 subtree2 = tree.get('child2')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	177 #TODO: Only simulate errors on portions of fragment that will become reads.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	178 add_pcr_errors(subtree1, '+', len(raw_frag_full), args.pcr_error, args.indel_rate, args.ext_rate)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	179 add_pcr_errors(subtree2, '-', len(raw_frag_full), args.pcr_error, args.indel_rate, args.ext_rate)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	180 apply_pcr_errors(tree, raw_frag_full)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	181 fragments = get_final_fragments(tree)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	182 add_mutation_lists(tree, fragments, [])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	183
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	184 # Step 4: Introduce sequencing errors.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	185 for fragment in fragments.values():
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	186 for mutation in generate_mutations(args.read_len, args.seq_error, args.indel_rate,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	187 args.ext_rate):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	188 fragment['mutations'].append(mutation)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	189 fragment['seq'] = apply_mutation(mutation, fragment['seq'])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	190
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	191 # Print barcodes to log file.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	192 if args.barcodes:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	193 args.barcodes.write('{}-{}\t{}\t{}\n'.format(chrom, id_num, barcode1, barcode2_rc))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	194 # Print family.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	195 for frag_id in sorted(fragments.keys()):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	196 fragment = fragments[frag_id]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	197 read_id = '{}-{}-{}'.format(chrom, id_num, frag_id)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	198 # Print mutations to log file.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	199 if args.mutations:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	200 read1_muts = get_mutations_subset(fragment['mutations'], 0, args.read_len)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	201 read2_muts = get_mutations_subset(fragment['mutations'], 0, args.read_len, revcomp=True,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	202 seqlen=len(fragment['seq']))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	203 if fragment['strand'] == '-':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	204 read1_muts, read2_muts = read2_muts, read1_muts
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	205 log_mutations(args.mutations, read1_muts, read_id+'/1', chrom, start, stop)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	206 log_mutations(args.mutations, read2_muts, read_id+'/2', chrom, start, stop)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	207 frag_seq = fragment['seq']
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	208 read1_seq = frag_seq[:args.read_len]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	209 read2_seq = get_revcomp(frag_seq[len(frag_seq)-args.read_len:])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	210 if fragment['strand'] == '-':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	211 read1_seq, read2_seq = read2_seq, read1_seq
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	212 if args.out_format == 'fasta':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	213 out1.write('>{}\n{}\n'.format(read_id, read1_seq))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	214 out2.write('>{}\n{}\n'.format(read_id, read2_seq))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	215 elif args.out_format == 'fastq':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	216 out1.write('@{}\n{}\n+\n{}\n'.format(read_id, read1_seq, qual_line))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	217 out2.write('@{}\n{}\n+\n{}\n'.format(read_id, read2_seq, qual_line))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	218
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	219 finally:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	220 try:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	221 os.remove(tmpfile.name)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	222 except OSError:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	223 pass
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	224
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	225
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	226 def run_command(command, *kwargs):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	227 """Run a command and return the exit code.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	228 run_command('echo', 'hello')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	229 Will print the command to stderr before running, unless "silent" is set to True."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	230 command_strs = map(str, command)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	231 if not kwargs.get('silent'):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	232 sys.stderr.write('$ '+' '.join(command_strs)+'\n')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	233 devnull = open(os.devnull, 'w')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	234 try:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	235 exit_status = subprocess.call(map(str, command), stderr=devnull)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	236 except OSError:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	237 exit_status = None
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	238 finally:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	239 devnull.close()
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	240 return exit_status
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	241
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	242
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	243 def extend_dist(raw_dist, exponent=1.25, min_prob=0.00001, max_len_mult=2):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	244 """Add an exponentially decreasing tail to the distribution.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	245 It takes the final value in the distribution and keeps dividing it by
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	246 "exponent", adding each new value to the end. It will not add probabilities
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	247 smaller than "min_prob" or extend the length of the list by more than
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	248 "max_len_mult" times."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	249 extended_dist = list(raw_dist)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	250 final_sum = sum(raw_dist)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	251 value = raw_dist[-1]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	252 value /= exponent
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	253 while value/final_sum >= min_prob and len(extended_dist) < len(raw_dist)*max_len_mult:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	254 extended_dist.append(value)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	255 final_sum += value
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	256 value /= exponent
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	257 return extended_dist
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	258
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	259
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	260 def compile_dist(raw_dist):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	261 """Turn the human-readable list of probabilities defined at the top into
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	262 proportional probabilities.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	263 E.g. [10, 5, 5] -> [0.5, 0.75, 1.0]"""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	264 proportional_dist = []
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	265 final_sum = sum(raw_dist)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	266 current_sum = 0
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	267 for magnitude in raw_dist:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	268 current_sum += magnitude
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	269 proportional_dist.append(current_sum/final_sum)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	270 return proportional_dist
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	271
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	272
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	273 def parse_read_id(read_id):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	274 match = re.search(WGSIM_ID_REGEX, read_id)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	275 if match:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	276 chrom = match.group(1)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	277 start = match.group(2)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	278 stop = match.group(3)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	279 id_num = match.group(4)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	280 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	281 chrom, id_num, start, stop = read_id, None, None, None
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	282 return chrom, id_num, start, stop
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	283
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	284
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	285 #TODO: Clean up "mutation" vs "error" terminology.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	286 def generate_mutations(seq_len, error_rate, indel_rate, extension_rate):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	287 """Generate all the mutations that occur over the length of a sequence."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	288 i = 0
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	289 while i <= seq_len:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	290 if random.random() < error_rate:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	291 mtype, alt = make_mutation(indel_rate, extension_rate)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	292 # Allow mutation after the last base only if it's an insertion.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	293 if i < seq_len or mtype == 'ins':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	294 yield {'coord':i, 'type':mtype, 'alt':alt}
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	295 # Compensate for length variations to keep i tracking the original read's base coordinates.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	296 if mtype == 'ins':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	297 i += len(alt)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	298 elif mtype == 'del':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	299 i -= alt
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	300 i += 1
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	301
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	302
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	303 def make_mutation(indel_rate, extension_rate):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	304 """Simulate a random mutation."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	305 # Is it an indel?
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	306 rand = random.random()
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	307 if rand < indel_rate:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	308 # Is it an insertion or deletion? Decide, then initialize it.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	309 # Re-use the random number from above. Just check if it's in the lower or upper half of the
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	310 # range from 0 to indel_rate.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	311 if rand < indel_rate/2:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	312 mtype = 'del'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	313 alt = 1
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	314 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	315 mtype = 'ins'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	316 alt = get_rand_base()
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	317 # Extend the indel as long as the extension rate allows.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	318 while random.random() < extension_rate:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	319 if mtype == 'ins':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	320 alt += get_rand_base()
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	321 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	322 alt += 1
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	323 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	324 # What is the new base for the SNV?
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	325 mtype = 'snv'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	326 alt = get_rand_base()
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	327 return mtype, alt
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	328
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	329
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	330 def get_rand_base(bases='ACGT'):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	331 return random.choice(bases)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	332
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	333
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	334 def get_rand_seq(seq_len):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	335 return ''.join([get_rand_base() for i in range(seq_len)])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	336
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	337
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	338 def get_revcomp(seq):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	339 return seq.translate(REVCOMP_TABLE)[::-1]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	340
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	341
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	342 def apply_mutation(mut, seq):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	343 i = mut['coord']
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	344 if mut['type'] == 'snv':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	345 # Replace the base at "coord".
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	346 new_seq = seq[:i] + mut['alt'] + seq[i+1:]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	347 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	348 # Indels are handled by inserting or deleting bases starting before the base at "coord".
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	349 # This goes agains the VCF convention, but it allows deleting the first and last base, as well
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	350 # as inserting before and after the sequence without as much special-casing.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	351 if mut['type'] == 'ins':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	352 # Example: 'ACGTACGT' + ins 'GC' at 4 = 'ACGTGCACGT'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	353 new_seq = seq[:i] + mut['alt'] + seq[i:]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	354 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	355 # Example: 'ACGTACGT' + del 2 at 4 = 'ACGTGT'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	356 new_seq = seq[:i] + seq[i+mut['alt']:]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	357 return new_seq
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	358
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	359
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	360 def get_mutations_subset(mutations_old, start, length, revcomp=False, seqlen=None):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	361 """Get a list of the input mutations which are within a certain region.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	362 The output list maintains the order in the input list, only filtering out
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	363 mutations outside the specified region.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	364 "start" is the start of the region (0-based). If revcomp, this start should be
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	365 in the coordinate system of the reverse-complemented sequence.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	366 "length" is the length of the region.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	367 "revcomp" causes the mutations to be converted to their reverse complements, and
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	368 the "start" to refer to the reverse complement sequence's coordinates. The order
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	369 of the mutations is unchanged, though.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	370 "seqlen" is the length of the sequence the mutations occurred in. This is only
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	371 needed when revcomp is True, to convert coordinates to the reverse complement
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	372 coordinate system."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	373 stop = start + length
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	374 mutations_new = []
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	375 for mutation in mutations_old:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	376 if revcomp:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	377 mutation = get_mutation_revcomp(mutation, seqlen)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	378 if start <= mutation['coord'] < stop:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	379 mutations_new.append(mutation)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	380 elif mutation['coord'] == stop and mutation['type'] == 'ins':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	381 # Allow insertions at the last coordinate.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	382 mutations_new.append(mutation)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	383 return mutations_new
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	384
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	385
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	386 def get_mutation_revcomp(mut, seqlen):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	387 """Convert a mutation to its reverse complement.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	388 "seqlen" is the length of the sequence the mutation is being applied to. Needed
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	389 to convert the coordinate to a coordinate system starting at the end of the
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	390 sequence."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	391 mut_rc = {'type':mut['type']}
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	392 if mut['type'] == 'snv':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	393 mut_rc['coord'] = seqlen - mut['coord'] - 1
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	394 mut_rc['alt'] = get_revcomp(mut['alt'])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	395 elif mut['type'] == 'ins':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	396 mut_rc['coord'] = seqlen - mut['coord']
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	397 mut_rc['alt'] = get_revcomp(mut['alt'])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	398 elif mut['type'] == 'del':
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	399 mut_rc['coord'] = seqlen - mut['coord'] - mut['alt']
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	400 mut_rc['alt'] = mut['alt']
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	401 return mut_rc
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	402
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	403
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	404 def log_mutations(mutfile, mutations, read_id, chrom, start, stop):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	405 for mutation in mutations:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	406 mutfile.write('{read_id}\t{chrom}\t{start}\t{stop}\t{coord}\t{type}\t{alt}\n'
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	407 .format(read_id=read_id, chrom=chrom, start=start, stop=stop, **mutation))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	408
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	409
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	410 def add_pcr_errors(subtree, strand, read_len, error_rate, indel_rate, extension_rate):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	411 """Add simulated PCR errors to a node in a tree and all its descendants."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	412 # Note: The errors are intended as "errors made in creating this molecule", so don't apply this to
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	413 # the root node, since that is supposed to be the original, unaltered molecule.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	414 # Go down the subtree and simulate errors in creating each fragment.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	415 # Process all the first-child descendants of the original node in a loop, and recursively call
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	416 # this function to process all second children.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	417 node = subtree
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	418 while node:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	419 node['strand'] = strand
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	420 node['errors'] = list(generate_mutations(read_len, error_rate, indel_rate, extension_rate))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	421 add_pcr_errors(node.get('child2'), strand, read_len, error_rate, indel_rate, extension_rate)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	422 node = node.get('child1')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	423
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	424
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	425 def apply_pcr_errors(subtree, seq):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	426 node = subtree
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	427 while node:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	428 for error in node.get('errors', ()):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	429 seq = apply_mutation(error, seq)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	430 if 'child1' not in node:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	431 node['seq'] = seq
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	432 apply_pcr_errors(node.get('child2'), seq)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	433 node = node.get('child1')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	434
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	435
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	436 def get_final_fragments(tree):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	437 """Walk to the leaf nodes of the tree and get the post-PCR sequences of all the fragments.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	438 Returns a dict mapping fragment id number to a dict representing the fragment. Its only two keys
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	439 are 'seq' (the final sequence) and 'strand' ('+' or '-')."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	440 fragments = {}
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	441 nodes = [tree]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	442 while nodes:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	443 node = nodes.pop()
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	444 child1 = node.get('child1')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	445 if child1:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	446 nodes.append(child1)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	447 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	448 fragments[node['id']] = {'seq':node['seq'], 'strand':node['strand']}
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	449 child2 = node.get('child2')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	450 if child2:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	451 nodes.append(child2)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	452 return fragments
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	453
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	454
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	455 def add_mutation_lists(subtree, fragments, mut_list1):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	456 """Compile the list of mutations that each fragment has undergone in PCR.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	457 To call from the root, give [] as "mut_list1" and a dict mapping all existing node id's to a dict
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	458 as "fragments". Instead of returning the data, this will add a 'mutations' key to the dict for
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	459 each fragment, mapping it to a list of PCR mutations that occurred in the lineage of the fragment,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	460 in chronological order."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	461 node = subtree
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	462 while node:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	463 mut_list1.extend(node.get('errors', ()))
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	464 if 'child1' not in node:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	465 fragments[node['id']]['mutations'] = mut_list1
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	466 if 'child2' in node:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	467 mut_list2 = copy.deepcopy(mut_list1)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	468 add_mutation_lists(node.get('child2'), fragments, mut_list2)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	469 node = node.get('child1')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	470
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	471
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	472 def check_tree_balance(subtree):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	473 """Find all points in the tree where the cycles of sibling nodes is unequal, and
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	474 return the maximum difference."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	475 node = subtree
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	476 if node:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	477 child1 = node.get('child1')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	478 child2 = node.get('child2')
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	479 if child1 and child2:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	480 diff = abs(child1['cycle'] - child2['cycle'])
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	481 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	482 diff = 0
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	483 diff_child1 = check_tree_balance(child1)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	484 diff_child2 = check_tree_balance(child2)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	485 return max(diff, diff_child1, diff_child2)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	486 else:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	487 return 0
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	488
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	489
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	490 def get_good_pcr_tree(n_reads, n_cycles, max_tries, max_diff=1):
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	491 """Return a single, balanced PCR tree from build_pcr_tree(), or fail if one cannot
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	492 be found in max_tries.
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	493 Compensate for bugs in build_pcr_tree() that sometimes result in multiple trees,
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	494 or trees with siblings from different cycles."""
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	495 tries = 0
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	496 while tries <= max_tries:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	497 trees = build_pcr_tree(n_reads, n_cycles)
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	498 if len(trees) == 1 and check_tree_balance(trees[0]) <= max_diff:
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	499 return trees[0]
e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty nick parents: diff changeset	500 tries += 1

18

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

1 #!/usr/bin/env python

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

2 from __future__ import division

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

3 from __future__ import print_function

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

4 import re

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

5 import os

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

6 import sys

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

7 import copy

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

8 import numpy

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

9 import bisect

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

10 import random

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

11 import string

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

12 import numbers

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

13 import tempfile

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

14 import argparse

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

15 import subprocess

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

16 import fastqreader

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

17

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

18 REVCOMP_TABLE = string.maketrans('acgtrymkbdhvACGTRYMKBDHV', 'tgcayrkmvhdbTGCAYRKMVHDB')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

19 WGSIM_ID_REGEX = r'^(.+)_(\d+)_(\d+)_\d+:\d+:\d+_\d+:\d+:\d+_([0-9a-f]+)/[12]$'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

20 ARG_DEFAULTS = {'read_len':100, 'frag_len':400, 'n_frags':1000, 'out_format':'fasta',

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

21 'seq_error':0.001, 'pcr_error':0.001, 'cycles':25, 'indel_rate':0.15,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

22 'ext_rate':0.3, 'seed':None, 'invariant':'TGACT', 'bar_len':12, 'fastq_qual':'I'}

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

23 USAGE = "%(prog)s [options]"

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

24 DESCRIPTION = """Simulate a duplex sequencing experiment."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

25

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

26 RAW_DISTRIBUTION = (

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

27 # 0 1 2 3 4 5 6 7 8 9

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

28 # Low singletons, but then constant drop-off. From pML113 (see 2015-09-28 report).

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

29 # 0, 100, 36, 31, 27, 22, 17, 12, 7, 4.3,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

30 #2.4, 1.2, 0.6, 0.3, 0.2, 0.15, 0.1, 0.07, 0.05, 0.03,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

31 # High singletons, but then a second peak around 10. From Christine plasmid (2015-10-06 report).

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

32 # 0, 100, 5.24, 3.67, 3.50, 3.67, 3.85, 4.02, 4.11, 4.20,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

33 # 4.17, 4.10, 4.00, 3.85, 3.69, 3.55, 3.38, 3.15, 2.92, 2.62,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

34 # 2.27, 2.01, 1.74, 1.56, 1.38, 1.20, 1.02, 0.85,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

35 # Same as above, but low singletons, 2's, and 3's (rely on errors to fill out those).

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

36 0, 1, 2, 3, 3.50, 3.67, 3.85, 4.02, 4.11, 4.20,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

37 4.17, 4.10, 4.00, 3.85, 3.69, 3.55, 3.38, 3.15, 2.92, 2.62,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

38 2.27, 2.01, 1.74, 1.56, 1.38, 1.20, 1.02, 0.85,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

39 )

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

40

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

41

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

42 def main(argv):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

43

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

44 parser = argparse.ArgumentParser(description=DESCRIPTION)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

45 parser.set_defaults(**ARG_DEFAULTS)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

46

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

47 parser.add_argument('ref', metavar='ref.fa', nargs='?',

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

48 help='Reference sequence. Omit if giving --frag-file.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

49 parser.add_argument('out1', type=argparse.FileType('w'),

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

50 help='Write final mate 1 reads to this file.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

51 parser.add_argument('out2', type=argparse.FileType('w'),

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

52 help='Write final mate 2 reads to this file.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

53 parser.add_argument('-o', '--out-format', choices=('fastq', 'fasta'))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

54 parser.add_argument('--stdout', action='store_true',

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

55 help='Print interleaved output reads to stdout.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

56 parser.add_argument('-m', '--mutations', type=argparse.FileType('w'),

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

57 help='Write a log of the PCR and sequencing errors introduced to this file. Will overwrite any '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

58 'existing file at this path.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

59 parser.add_argument('-b', '--barcodes', type=argparse.FileType('w'),

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

60 help='Write a log of which barcodes were ligated to which fragments. Will overwrite any '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

61 'existing file at this path.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

62 parser.add_argument('--frag-file',

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

63 help='The path of the FASTQ file of fragments. If --ref is given, these will be generated with '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

64 'wgsim and kept (normally a temporary file is used, then deleted). Note: the file will be '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

65 'overwritten! If --ref is not given, then this should be a file of already generated '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

66 'fragments, and they will be used instead of generating new ones.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

67 parser.add_argument('-Q', '--fastq-qual',

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

68 help='The quality score to assign to all bases in FASTQ output. Give a character or PHRED '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

69 'score (integer). A PHRED score will be converted using the Sanger offset (33). Default: '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

70 '"%(default)s"')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

71 parser.add_argument('-S', '--seed', type=int,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

72 help='Random number generator seed. By default, a random, 32-bit seed will be generated and '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

73 'logged to stdout.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

74 params = parser.add_argument_group('simulation parameters')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

75 params.add_argument('-n', '--n-frags', type=int,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

76 help='The number of original fragment molecules to simulate. The final number of reads will be '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

77 'this multiplied by the average number of reads per family. If you provide fragments with '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

78 '--frag-file, the script will still only read in the number specified here. Default: '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

79 '%(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

80 params.add_argument('-r', '--read-len', type=int,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

81 help='Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

82 params.add_argument('-f', '--frag-len', type=int,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

83 help='Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

84 params.add_argument('-s', '--seq-error', type=float,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

85 help='Sequencing error rate per base (0-1 proportion, not percent). Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

86 params.add_argument('-p', '--pcr-error', type=float,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

87 help='PCR error rate per base (0-1 proportion, not percent). Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

88 params.add_argument('-c', '--cycles', type=int,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

89 help='Number of PCR cycles to simulate. Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

90 params.add_argument('-i', '--indel-rate', type=float,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

91 help='Fraction of errors which are indels. Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

92 params.add_argument('-E', '--extension-rate', dest='ext_rate', type=float,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

93 help='Probability an indel is extended. Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

94 params.add_argument('-B', '--bar-len', type=int,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

95 help='Length of the barcodes to generate. Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

96 params.add_argument('-I', '--invariant',

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

97 help='The invariant linker sequence between the barcode and sample sequence in each read. '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

98 'Default: %(default)s')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

99

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

100 # Parse and interpret arguments.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

101 args = parser.parse_args(argv[1:])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

102 assert args.ref or args.frag_file, 'You must provide either a reference or fragments file.'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

103 if args.seed is None:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

104 seed = random.randint(0, 2**31-1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

105 sys.stderr.write('seed: {}\n'.format(seed))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

106 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

107 seed = args.seed

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

108 random.seed(seed)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

109 if args.stdout:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

110 out1 = sys.stdout

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

111 out2 = sys.stdout

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

112 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

113 out1 = args.out1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

114 out2 = args.out2

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

115 if isinstance(args.fastq_qual, numbers.Integral):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

116 assert args.fastq_qual >= 0, '--fastq-qual cannot be negative.'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

117 fastq_qual = chr(args.fastq_qual + 33)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

118 elif isinstance(args.fastq_qual, basestring):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

119 assert len(args.fastq_qual) == 1, '--fastq-qual cannot be more than a single character.'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

120 fastq_qual = args.fastq_qual

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

121 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

122 raise AssertionError('--fastq-qual must be a positive integer or single character.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

123 qual_line = fastq_qual * args.read_len

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

124

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

125 invariant_rc = get_revcomp(args.invariant)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

126

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

127 # Create a temporary director to do our work in. Then work inside a try so we can finally remove

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

128 # the directory no matter what exceptions are encountered.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

129 tmpfile = tempfile.NamedTemporaryFile(prefix='wgdsim.frags.')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

130 tmpfile.close()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

131 try:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

132 # Step 1: Use wgsim to create fragments from the reference.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

133 if args.frag_file:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

134 frag_file = args.frag_file

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

135 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

136 frag_file = tmpfile.name

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

137 if args.ref and os.path.isfile(args.ref) and os.path.getsize(args.ref):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

138 #TODO: Check exit status

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

139 #TODO: Check for wgsim on the PATH.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

140 # Set error and mutation rates to 0 to just slice sequences out of the reference without

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

141 # modification.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

142 run_command('wgsim', '-e', '0', '-r', '0', '-d', '0', '-R', args.indel_rate, '-S', seed,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

143 '-N', args.n_frags, '-X', args.ext_rate, '-1', args.frag_len,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

144 args.ref, frag_file, os.devnull)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

145

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

146 # NOTE: Coordinates here are 0-based (0 is the first base in the sequence).

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

147 extended_dist = extend_dist(RAW_DISTRIBUTION)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

148 proportional_dist = compile_dist(extended_dist)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

149 n_frags = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

150 for raw_fragment in fastqreader.FastqReadGenerator(frag_file):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

151 n_frags += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

152 if n_frags > args.n_frags:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

153 break

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

154 chrom, id_num, start, stop = parse_read_id(raw_fragment.id)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

155 barcode1 = get_rand_seq(args.bar_len)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

156 barcode2 = get_rand_seq(args.bar_len)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

157 barcode2_rc = get_revcomp(barcode2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

158 raw_frag_full = barcode1 + args.invariant + raw_fragment.seq + invariant_rc + barcode2

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

159

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

160 # Step 2: Determine how many reads to produce from each fragment.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

161 # - Use random.random() and divide the range 0-1 into segments of sizes proportional to

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

162 # the likelihood of each family size.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

163 # bisect.bisect() finds where an element belongs in a sorted list, returning the index.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

164 # proportional_dist is just such a sorted list, with values from 0 to 1.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

165 n_reads = bisect.bisect(proportional_dist, random.random())

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

166

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

167 # Step 3: Introduce PCR errors.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

168 # - Determine the mutations and their frequencies.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

169 # - Could get frequency from the cycle of PCR it occurs in.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

170 # - Important to have PCR errors shared between reads.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

171 # - For each read, determine which mutations it contains.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

172 # - Use random.random() < mut_freq.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

173 tree = get_good_pcr_tree(n_reads, args.cycles, 1000, max_diff=1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

174 # Add errors to all children of original fragment.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

175 subtree1 = tree.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

176 subtree2 = tree.get('child2')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

177 #TODO: Only simulate errors on portions of fragment that will become reads.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

178 add_pcr_errors(subtree1, '+', len(raw_frag_full), args.pcr_error, args.indel_rate, args.ext_rate)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

179 add_pcr_errors(subtree2, '-', len(raw_frag_full), args.pcr_error, args.indel_rate, args.ext_rate)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

180 apply_pcr_errors(tree, raw_frag_full)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

181 fragments = get_final_fragments(tree)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

182 add_mutation_lists(tree, fragments, [])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

183

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

184 # Step 4: Introduce sequencing errors.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

185 for fragment in fragments.values():

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

186 for mutation in generate_mutations(args.read_len, args.seq_error, args.indel_rate,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

187 args.ext_rate):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

188 fragment['mutations'].append(mutation)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

189 fragment['seq'] = apply_mutation(mutation, fragment['seq'])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

190

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

191 # Print barcodes to log file.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

192 if args.barcodes:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

193 args.barcodes.write('{}-{}\t{}\t{}\n'.format(chrom, id_num, barcode1, barcode2_rc))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

194 # Print family.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

195 for frag_id in sorted(fragments.keys()):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

196 fragment = fragments[frag_id]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

197 read_id = '{}-{}-{}'.format(chrom, id_num, frag_id)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

198 # Print mutations to log file.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

199 if args.mutations:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

200 read1_muts = get_mutations_subset(fragment['mutations'], 0, args.read_len)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

201 read2_muts = get_mutations_subset(fragment['mutations'], 0, args.read_len, revcomp=True,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

202 seqlen=len(fragment['seq']))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

203 if fragment['strand'] == '-':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

204 read1_muts, read2_muts = read2_muts, read1_muts

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

205 log_mutations(args.mutations, read1_muts, read_id+'/1', chrom, start, stop)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

206 log_mutations(args.mutations, read2_muts, read_id+'/2', chrom, start, stop)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

207 frag_seq = fragment['seq']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

208 read1_seq = frag_seq[:args.read_len]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

209 read2_seq = get_revcomp(frag_seq[len(frag_seq)-args.read_len:])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

210 if fragment['strand'] == '-':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

211 read1_seq, read2_seq = read2_seq, read1_seq

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

212 if args.out_format == 'fasta':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

213 out1.write('>{}\n{}\n'.format(read_id, read1_seq))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

214 out2.write('>{}\n{}\n'.format(read_id, read2_seq))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

215 elif args.out_format == 'fastq':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

216 out1.write('@{}\n{}\n+\n{}\n'.format(read_id, read1_seq, qual_line))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

217 out2.write('@{}\n{}\n+\n{}\n'.format(read_id, read2_seq, qual_line))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

218

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

219 finally:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

220 try:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

221 os.remove(tmpfile.name)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

222 except OSError:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

223 pass

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

224

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

225

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

226 def run_command(*command, **kwargs):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

227 """Run a command and return the exit code.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

228 run_command('echo', 'hello')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

229 Will print the command to stderr before running, unless "silent" is set to True."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

230 command_strs = map(str, command)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

231 if not kwargs.get('silent'):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

232 sys.stderr.write('$ '+' '.join(command_strs)+'\n')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

233 devnull = open(os.devnull, 'w')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

234 try:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

235 exit_status = subprocess.call(map(str, command), stderr=devnull)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

236 except OSError:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

237 exit_status = None

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

238 finally:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

239 devnull.close()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

240 return exit_status

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

241

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

242

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

243 def extend_dist(raw_dist, exponent=1.25, min_prob=0.00001, max_len_mult=2):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

244 """Add an exponentially decreasing tail to the distribution.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

245 It takes the final value in the distribution and keeps dividing it by

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

246 "exponent", adding each new value to the end. It will not add probabilities

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

247 smaller than "min_prob" or extend the length of the list by more than

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

248 "max_len_mult" times."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

249 extended_dist = list(raw_dist)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

250 final_sum = sum(raw_dist)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

251 value = raw_dist[-1]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

252 value /= exponent

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

253 while value/final_sum >= min_prob and len(extended_dist) < len(raw_dist)*max_len_mult:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

254 extended_dist.append(value)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

255 final_sum += value

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

256 value /= exponent

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

257 return extended_dist

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

258

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

259

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

260 def compile_dist(raw_dist):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

261 """Turn the human-readable list of probabilities defined at the top into

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

262 proportional probabilities.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

263 E.g. [10, 5, 5] -> [0.5, 0.75, 1.0]"""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

264 proportional_dist = []

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

265 final_sum = sum(raw_dist)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

266 current_sum = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

267 for magnitude in raw_dist:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

268 current_sum += magnitude

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

269 proportional_dist.append(current_sum/final_sum)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

270 return proportional_dist

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

271

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

272

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

273 def parse_read_id(read_id):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

274 match = re.search(WGSIM_ID_REGEX, read_id)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

275 if match:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

276 chrom = match.group(1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

277 start = match.group(2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

278 stop = match.group(3)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

279 id_num = match.group(4)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

280 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

281 chrom, id_num, start, stop = read_id, None, None, None

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

282 return chrom, id_num, start, stop

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

283

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

284

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

285 #TODO: Clean up "mutation" vs "error" terminology.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

286 def generate_mutations(seq_len, error_rate, indel_rate, extension_rate):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

287 """Generate all the mutations that occur over the length of a sequence."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

288 i = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

289 while i <= seq_len:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

290 if random.random() < error_rate:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

291 mtype, alt = make_mutation(indel_rate, extension_rate)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

292 # Allow mutation after the last base only if it's an insertion.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

293 if i < seq_len or mtype == 'ins':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

294 yield {'coord':i, 'type':mtype, 'alt':alt}

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

295 # Compensate for length variations to keep i tracking the original read's base coordinates.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

296 if mtype == 'ins':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

297 i += len(alt)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

298 elif mtype == 'del':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

299 i -= alt

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

300 i += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

301

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

302

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

303 def make_mutation(indel_rate, extension_rate):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

304 """Simulate a random mutation."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

305 # Is it an indel?

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

306 rand = random.random()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

307 if rand < indel_rate:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

308 # Is it an insertion or deletion? Decide, then initialize it.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

309 # Re-use the random number from above. Just check if it's in the lower or upper half of the

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

310 # range from 0 to indel_rate.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

311 if rand < indel_rate/2:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

312 mtype = 'del'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

313 alt = 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

314 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

315 mtype = 'ins'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

316 alt = get_rand_base()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

317 # Extend the indel as long as the extension rate allows.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

318 while random.random() < extension_rate:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

319 if mtype == 'ins':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

320 alt += get_rand_base()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

321 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

322 alt += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

323 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

324 # What is the new base for the SNV?

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

325 mtype = 'snv'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

326 alt = get_rand_base()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

327 return mtype, alt

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

328

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

329

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

330 def get_rand_base(bases='ACGT'):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

331 return random.choice(bases)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

332

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

333

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

334 def get_rand_seq(seq_len):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

335 return ''.join([get_rand_base() for i in range(seq_len)])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

336

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

337

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

338 def get_revcomp(seq):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

339 return seq.translate(REVCOMP_TABLE)[::-1]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

340

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

341

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

342 def apply_mutation(mut, seq):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

343 i = mut['coord']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

344 if mut['type'] == 'snv':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

345 # Replace the base at "coord".

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

346 new_seq = seq[:i] + mut['alt'] + seq[i+1:]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

347 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

348 # Indels are handled by inserting or deleting bases starting *before* the base at "coord".

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

349 # This goes agains the VCF convention, but it allows deleting the first and last base, as well

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

350 # as inserting before and after the sequence without as much special-casing.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

351 if mut['type'] == 'ins':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

352 # Example: 'ACGTACGT' + ins 'GC' at 4 = 'ACGTGCACGT'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

353 new_seq = seq[:i] + mut['alt'] + seq[i:]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

354 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

355 # Example: 'ACGTACGT' + del 2 at 4 = 'ACGTGT'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

356 new_seq = seq[:i] + seq[i+mut['alt']:]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

357 return new_seq

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

358

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

359

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

360 def get_mutations_subset(mutations_old, start, length, revcomp=False, seqlen=None):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

361 """Get a list of the input mutations which are within a certain region.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

362 The output list maintains the order in the input list, only filtering out

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

363 mutations outside the specified region.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

364 "start" is the start of the region (0-based). If revcomp, this start should be

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

365 in the coordinate system of the reverse-complemented sequence.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

366 "length" is the length of the region.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

367 "revcomp" causes the mutations to be converted to their reverse complements, and

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

368 the "start" to refer to the reverse complement sequence's coordinates. The order

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

369 of the mutations is unchanged, though.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

370 "seqlen" is the length of the sequence the mutations occurred in. This is only

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

371 needed when revcomp is True, to convert coordinates to the reverse complement

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

372 coordinate system."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

373 stop = start + length

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

374 mutations_new = []

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

375 for mutation in mutations_old:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

376 if revcomp:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

377 mutation = get_mutation_revcomp(mutation, seqlen)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

378 if start <= mutation['coord'] < stop:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

379 mutations_new.append(mutation)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

380 elif mutation['coord'] == stop and mutation['type'] == 'ins':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

381 # Allow insertions at the last coordinate.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

382 mutations_new.append(mutation)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

383 return mutations_new

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

384

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

385

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

386 def get_mutation_revcomp(mut, seqlen):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

387 """Convert a mutation to its reverse complement.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

388 "seqlen" is the length of the sequence the mutation is being applied to. Needed

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

389 to convert the coordinate to a coordinate system starting at the end of the

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

390 sequence."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

391 mut_rc = {'type':mut['type']}

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

392 if mut['type'] == 'snv':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

393 mut_rc['coord'] = seqlen - mut['coord'] - 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

394 mut_rc['alt'] = get_revcomp(mut['alt'])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

395 elif mut['type'] == 'ins':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

396 mut_rc['coord'] = seqlen - mut['coord']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

397 mut_rc['alt'] = get_revcomp(mut['alt'])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

398 elif mut['type'] == 'del':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

399 mut_rc['coord'] = seqlen - mut['coord'] - mut['alt']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

400 mut_rc['alt'] = mut['alt']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

401 return mut_rc

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

402

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

403

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

404 def log_mutations(mutfile, mutations, read_id, chrom, start, stop):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

405 for mutation in mutations:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

406 mutfile.write('{read_id}\t{chrom}\t{start}\t{stop}\t{coord}\t{type}\t{alt}\n'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

407 .format(read_id=read_id, chrom=chrom, start=start, stop=stop, **mutation))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

408

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

409

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

410 def add_pcr_errors(subtree, strand, read_len, error_rate, indel_rate, extension_rate):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

411 """Add simulated PCR errors to a node in a tree and all its descendants."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

412 # Note: The errors are intended as "errors made in creating this molecule", so don't apply this to

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

413 # the root node, since that is supposed to be the original, unaltered molecule.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

414 # Go down the subtree and simulate errors in creating each fragment.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

415 # Process all the first-child descendants of the original node in a loop, and recursively call

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

416 # this function to process all second children.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

417 node = subtree

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

418 while node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

419 node['strand'] = strand

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

420 node['errors'] = list(generate_mutations(read_len, error_rate, indel_rate, extension_rate))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

421 add_pcr_errors(node.get('child2'), strand, read_len, error_rate, indel_rate, extension_rate)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

422 node = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

423

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

424

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

425 def apply_pcr_errors(subtree, seq):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

426 node = subtree

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

427 while node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

428 for error in node.get('errors', ()):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

429 seq = apply_mutation(error, seq)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

430 if 'child1' not in node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

431 node['seq'] = seq

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

432 apply_pcr_errors(node.get('child2'), seq)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

433 node = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

434

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

435

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

436 def get_final_fragments(tree):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

437 """Walk to the leaf nodes of the tree and get the post-PCR sequences of all the fragments.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

438 Returns a dict mapping fragment id number to a dict representing the fragment. Its only two keys

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

439 are 'seq' (the final sequence) and 'strand' ('+' or '-')."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

440 fragments = {}

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

441 nodes = [tree]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

442 while nodes:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

443 node = nodes.pop()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

444 child1 = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

445 if child1:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

446 nodes.append(child1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

447 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

448 fragments[node['id']] = {'seq':node['seq'], 'strand':node['strand']}

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

449 child2 = node.get('child2')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

450 if child2:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

451 nodes.append(child2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

452 return fragments

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

453

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

454

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

455 def add_mutation_lists(subtree, fragments, mut_list1):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

456 """Compile the list of mutations that each fragment has undergone in PCR.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

457 To call from the root, give [] as "mut_list1" and a dict mapping all existing node id's to a dict

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

458 as "fragments". Instead of returning the data, this will add a 'mutations' key to the dict for

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

459 each fragment, mapping it to a list of PCR mutations that occurred in the lineage of the fragment,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

460 in chronological order."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

461 node = subtree

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

462 while node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

463 mut_list1.extend(node.get('errors', ()))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

464 if 'child1' not in node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

465 fragments[node['id']]['mutations'] = mut_list1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

466 if 'child2' in node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

467 mut_list2 = copy.deepcopy(mut_list1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

468 add_mutation_lists(node.get('child2'), fragments, mut_list2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

469 node = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

470

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

471

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

472 def check_tree_balance(subtree):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

473 """Find all points in the tree where the cycles of sibling nodes is unequal, and

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

474 return the maximum difference."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

475 node = subtree

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

476 if node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

477 child1 = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

478 child2 = node.get('child2')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

479 if child1 and child2:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

480 diff = abs(child1['cycle'] - child2['cycle'])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

481 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

482 diff = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

483 diff_child1 = check_tree_balance(child1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

484 diff_child2 = check_tree_balance(child2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

485 return max(diff, diff_child1, diff_child2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

486 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

487 return 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

488

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

489

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

490 def get_good_pcr_tree(n_reads, n_cycles, max_tries, max_diff=1):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

491 """Return a single, balanced PCR tree from build_pcr_tree(), or fail if one cannot

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

492 be found in max_tries.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

493 Compensate for bugs in build_pcr_tree() that sometimes result in multiple trees,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

494 or trees with siblings from different cycles."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

495 tries = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

496 while tries <= max_tries:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

497 trees = build_pcr_tree(n_reads, n_cycles)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

498 if len(trees) == 1 and check_tree_balance(trees[0]) <= max_diff:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

499 return trees[0]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

500 tries += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

501 raise AssertionError('Could not generate a single, balanced tree! (tried {} times)'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

502 .format(max_tries))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

503

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

504

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

505 def build_pcr_tree(n_reads, n_cycles):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

506 """Create a simulated descent lineage of how all the final PCR fragments are related.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

507 Each node represents a fragment molecule at one stage of PCR. Each node is a dict containing the

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

508 fragment's children (other nodes) ('child1' and 'child2'), the PCR cycle number ('cycle'), and,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

509 at the leaves, a unique id number for each final fragment.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

510 Returns a list of root nodes. Usually there will only be one, but about 1-3% of the time it fails

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

511 to unify the subtrees and results in a broken tree.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

512 """

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

513 #TODO: Make it always return a single tree.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

514 # Begin a branch for each of the fragments. These are the leaf nodes. We'll work backward from

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

515 # these, simulating the points at which they share ancestors, eventually coalescing into the

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

516 # single (root) ancestor.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

517 branches = []

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

518 for frag_id in range(n_reads):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

519 branches.append({'cycle':n_cycles-1, 'id':frag_id})

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

520 # Build up all the branches in parallel. Start from the second-to-last PCR cycle.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

521 for cycle in reversed(range(n_cycles-1)):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

522 # Probability of 2 fragments sharing an ancestor at cycle c is 1/2^c.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

523 prob = 1/2**cycle

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

524 frag_i = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

525 while frag_i < len(branches):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

526 current_root = branches[frag_i]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

527 # Does the current fragment share this ancestor with any of the other fragments?

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

528 # numpy.random.binomial() is a fast way to simulate going through every other fragment and

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

529 # asking if random.random() < prob.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

530 shared = numpy.random.binomial(len(branches)-1, prob)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

531 if shared == 0:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

532 # No branch point here. Just add another level to the lineage.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

533 branches[frag_i] = {'cycle':cycle, 'child1':current_root}

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

534 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

535 # Pick a random other fragment to share this ancestor with.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

536 # Make a list of candidates to pick from.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

537 candidates = []

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

538 for candidate_i, candidate in enumerate(branches):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

539 # Don't include ourselves.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

540 if candidate is current_root:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

541 continue

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

542 # If it's at a cycle above us and it already has a child, skip it.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

543 if candidate['cycle'] == cycle and candidate.get('child2'):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

544 continue

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

545 candidates.append(candidate_i)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

546 if candidates:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

547 relative_i = random.choice(candidates)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

548 relative = branches[relative_i]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

549 # Have we already passed this fragmentfragment on this cycle?

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

550 if relative['cycle'] == cycle:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

551 # If we've already passed it, we're looking at the fragment's parent. We want the child.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

552 relative = relative['child1']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

553 # Join the lineages of our current fragment and the relative to a new parent.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

554 #TODO: Sometimes, we end up matching up subtrees of different depths. But the discrepancy

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

555 # is rarely greater than 1. Figure out why.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

556 # assert abs(current_root['cycle'] - relative['cycle']) < 3, ('cycle: {}, current_root: {},'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

557 # ' relative: {}, frag_i: {}, relative_i: {}, branches: {}, candidates: {}, shared: {}'

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

558 # .format(cycle, current_root['cycle'], relative['cycle'], frag_i, relative_i,

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

559 # len(branches), len(candidates), shared))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

560 branches[frag_i] = {'cycle':cycle, 'child1':current_root, 'child2':relative}

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

561 # Remove the relative from the list of lineages.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

562 del(branches[relative_i])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

563 if relative_i < frag_i:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

564 frag_i -= 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

565 frag_i += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

566 return branches

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

567

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

568

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

569 def get_depth(tree):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

570 depth = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

571 node = tree

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

572 while node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

573 depth += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

574 node = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

575 return depth

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

576

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

577

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

578 def convert_tree(tree_orig):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

579 # Let's operate on a copy only.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

580 tree = copy.deepcopy(tree_orig)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

581 # Turn the tree vertical.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

582 tree['line'] = 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

583 tree['children'] = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

584 levels = [[tree]]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

585 done = False

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

586 while not done:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

587 last_level = levels[-1]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

588 this_level = []

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

589 done = True

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

590 for node in last_level:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

591 for child_name in ('child1', 'child2'):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

592 child = node.get(child_name)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

593 if child:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

594 done = False

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

595 child['parent'] = node

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

596 child['branch'] = child['parent']['branch']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

597 if child_name == 'child2':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

598 child['branch'] += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

599 this_level.append(child)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

600 this_level.sort(key=lambda node: node['branch'])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

601 levels.append(this_level)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

602 return levels

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

603

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

604

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

605 def print_levels(levels):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

606 last_level = []

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

607 for level in levels:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

608 for node in level:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

609 child = 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

610 for parent in last_level:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

611 if parent.get('child2') is node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

612 child = 2

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

613 if child == 1:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

614 sys.stdout.write('| ')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

615 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

616 sys.stdout.write('\ ')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

617 last_level = level

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

618 print()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

619

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

620

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

621 def label_branches(tree):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

622 """Label each vertical branch (line of 'child1's) with an id number."""

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

623 counter = 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

624 tree['branch'] = counter

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

625 nodes = [tree]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

626 while nodes:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

627 node = nodes.pop(0)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

628 child1 = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

629 if child1:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

630 child1['branch'] = node['branch']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

631 nodes.append(child1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

632 child2 = node.get('child2')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

633 if child2:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

634 counter += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

635 child2['branch'] = counter

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

636 nodes.append(child2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

637

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

638

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

639 def print_tree(tree_orig):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

640 # We "write" strings to an output buffer instead of directly printing, so we can post-process the

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

641 # output. The buffer is a matrix of cells, each holding a string representing one element.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

642 lines = [[]]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

643 # Let's operate on a copy only.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

644 tree = copy.deepcopy(tree_orig)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

645 # Add some bookkeeping data.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

646 label_branches(tree)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

647 tree['level'] = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

648 branches = [tree]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

649 while branches:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

650 line = lines[-1]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

651 branch = branches.pop()

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

652 level = branch['level']

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

653 while level > 0:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

654 line.append(' ')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

655 level -= 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

656 node = branch

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

657 while node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

658 # Is it the root node? (Have we written anything yet?)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

659 if lines[0]:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

660 # Are we at the start of the line? (Is it only spaces so far?)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

661 if line[-1] == ' ':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

662 line.append('\-')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

663 elif line[-1].endswith('-'):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

664 line.append('=-')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

665 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

666 line.append('*-')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

667 child2 = node.get('child2')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

668 if child2:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

669 child2['level'] = node['level'] + 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

670 branches.append(child2)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

671 parent = node

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

672 node = node.get('child1')

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

673 if node:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

674 node['level'] = parent['level'] + 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

675 else:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

676 line.append(' {}'.format(parent['branch']))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

677 lines.append([])

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

678 # Post-process output: Add lines connecting branches to parents.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

679 x = 0

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

680 done = False

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

681 while not done:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

682 # Draw vertical lines upward from branch points.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

683 drawing = False

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

684 for line in reversed(lines):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

685 done = True

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

686 if x < len(line):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

687 done = False

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

688 cell = line[x]

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

689 if cell == '\-':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

690 drawing = True

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

691 elif cell == ' ' and drawing:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

692 line[x] = '| '

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

693 elif cell == '=-' and drawing:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

694 drawing = False

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

695 x += 1

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

696 # Print the final output.

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

697 for line in lines:

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

698 print(''.join(line))

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

699

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

700

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

701 def fail(message):

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

702 sys.stderr.write(message+"\n")

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

703 sys.exit(1)

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

704

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

705 if __name__ == '__main__':

e4d75f9efb90 planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty

nick

parents:

diff changeset

706 sys.exit(main(sys.argv))

Mercurial > repos > nick > duplex

annotate utils/sim.py @ 18:e4d75f9efb90 draft