<s> import sys <EOL> from collections import namedtuple <EOL> from seqio import fastaIterator <EOL> from itertools import dropwhile , count , repeat <EOL> import misc <EOL> import random <EOL> def insertion ( pos , seq ) : <EOL> return seq [ : pos ] + random . choice ( [ "<STR_LIT:A>" , "<STR_LIT:C>" , "<STR_LIT>" , "<STR_LIT:T>" ] ) + seq [ pos : ] <EOL> def deletion ( pos , seq ) : <EOL> return seq [ : pos ] + seq [ pos + <NUM_LIT:1> : ] <EOL> def mismatch ( pos , seq ) : <EOL> return seq [ : pos ] + random . choice ( [ "<STR_LIT:A>" , "<STR_LIT:C>" , "<STR_LIT>" , "<STR_LIT:T>" ] ) + seq [ pos + <NUM_LIT:1> : ] <EOL> if not len ( sys . argv ) == <NUM_LIT:5> : <EOL> print "<STR_LIT>" <EOL> sys . exit ( <NUM_LIT:1> ) <EOL> Chromosome = namedtuple ( "<STR_LIT>" , [ "<STR_LIT:name>" , "<STR_LIT>" ] ) <EOL> gfh = open ( sys . argv [ <NUM_LIT:1> ] ) <EOL> lfh = open ( sys . argv [ <NUM_LIT:2> ] ) <EOL> erate = float ( sys . argv [ <NUM_LIT:3> ] ) <EOL> rout = open ( sys . argv [ <NUM_LIT:4> ] + "<STR_LIT>" , "<STR_LIT:w>" ) <EOL> chromosomes = map ( lambda r : Chromosome . _make ( ( str ( r . name ) , str ( r . seq ) ) ) , fastaIterator ( gfh ) ) <EOL> chrom_lengths = map ( lambda c : len ( c . seq ) , chromosomes ) <EOL> genome_length = sum ( chrom_lengths ) <EOL> chrom_lengths_ivtf = map ( misc . accumulator ( <NUM_LIT:0> ) , map ( lambda x : float ( x ) / genome_length , chrom_lengths ) ) <EOL> count = <NUM_LIT:0> <EOL> for l in lfh : <EOL> readlen = int ( l . strip ( ) ) <EOL> while True : <EOL> U_c = random . random ( ) <EOL> chr_idx = misc . first_idx ( lambda x : x > U_c , chrom_lengths_ivtf ) <EOL> chromosome = chromosomes [ chr_idx ] <EOL> chromosome_len = chrom_lengths [ chr_idx ] <EOL> U_p = int ( random . random ( ) * chromosome_len ) <EOL> if ( U_p + readlen ) > ( chromosome_len - <NUM_LIT:1> ) : <EOL> continue <EOL> f = random . choice ( [ ( "<STR_LIT:f>" , lambda x : x ) , ( "<STR_LIT:r>" , misc . reverse_complement ) ] ) <EOL> sim_read = f [ <NUM_LIT:1> ] ( chromosome . seq [ U_p : U_p + readlen ] ) <EOL> sim_read_name = "<STR_LIT:>>" + chromosome . name + "<STR_LIT::>" + str ( U_p ) + "<STR_LIT:->" + str ( U_p + readlen ) + "<STR_LIT::>" + f [ <NUM_LIT:0> ] <EOL> num_errors = abs ( int ( random . normalvariate ( readlen * erate , readlen * erate * ( <NUM_LIT:1> - erate ) ) ) ) <EOL> for i in range ( num_errors ) : <EOL> position = int ( random . random ( ) * len ( sim_read ) ) <EOL> f = random . choice ( [ insertion , deletion , mismatch ] ) <EOL> sim_read = f ( position , sim_read ) <EOL> sim_read_name += "<STR_LIT>" % ( readlen , len ( sim_read ) , num_errors ) <EOL> rout . write ( "<STR_LIT>" % ( sim_read_name , sim_read ) ) <EOL> break <EOL> count += <NUM_LIT:1> <EOL> if count % <NUM_LIT> == <NUM_LIT:0> : <EOL> sys . stderr . write ( "<STR_LIT>" % count ) <EOL> gfh . close ( ) <EOL> lfh . close ( ) <EOL> rout . close ( ) </s>
