def find_23S(fastas, hmms, bit_thresh = float(<NUM_LIT:20>), length_thresh = <NUM_LIT>, masking = True, buffer = <NUM_LIT:0>):
<EOL>group2hmm = find_coordinates(hmms, bit_thresh)<EOL>for fasta in fastas:<EOL><INDENT>for seq in parse_fasta(fasta):<EOL><INDENT>id = seq[<NUM_LIT:0>].split('<STR_LIT:>>')[<NUM_LIT:1>].split()[<NUM_LIT:0>]<EOL>if id not in group2hmm:<EOL><INDENT>continue<EOL><DEDENT>seq[<NUM_LIT:1>] = seq[<NUM_LIT:1>].upper()<EOL>count = <NUM_LIT:0> <EOL>for group, info in list(group2hmm[id].items()):<EOL><INDENT>model, strand, coords, matches, gaps = info<EOL>ib = sum([i[<NUM_LIT:1>] - i[<NUM_LIT:0>] + <NUM_LIT:1> for i in gaps])<EOL>tl = coords[<NUM_LIT:1>] - coords[<NUM_LIT:0>] + <NUM_LIT:1><EOL>length = tl - ib<EOL>if length < length_thresh:<EOL><INDENT>continue <EOL><DEDENT>count += <NUM_LIT:1><EOL>ret_coords = [max([coords[<NUM_LIT:0>] - buffer, <NUM_LIT:1>]),min([coords[<NUM_LIT:1>] + buffer, len(seq[<NUM_LIT:1>])]), coords[<NUM_LIT:2>]]<EOL>buffer_ends = check_buffer(coords, len(seq[<NUM_LIT:1>]), buffer)<EOL>if masking is True:<EOL><INDENT>seq[<NUM_LIT:1>] = mask_sequence(seq[<NUM_LIT:1>], gaps)<EOL><DEDENT>S = seq[<NUM_LIT:1>][(ret_coords[<NUM_LIT:0>] - <NUM_LIT:1>):(ret_coords[<NUM_LIT:1>])]<EOL>inserts = [gap[<NUM_LIT:1>] - gap[<NUM_LIT:0>] + <NUM_LIT:1> for gap in gaps]<EOL>inserts.append('<STR_LIT:end>')<EOL>model_pos = '<STR_LIT:;>'.join(['<STR_LIT>' % (match[<NUM_LIT:2>], match[<NUM_LIT:3>], insert) for match, insert in zip(matches, inserts)])<EOL>header = '<STR_LIT>' %(seq[<NUM_LIT:0>], model, count, ret_coords[<NUM_LIT:0>], ret_coords[<NUM_LIT:1>], strand, tl, length, model_pos, buffer_ends[<NUM_LIT:0>], buffer_ends[<NUM_LIT:1>], ib)<EOL>if strand == '<STR_LIT:->':<EOL><INDENT>S = rc(['<STR_LIT>', S])[<NUM_LIT:1>]<EOL><DEDENT>yield [header, S]<EOL><DEDENT><DEDENT><DEDENT>
1) parse hmm output into dictionary (sequence must pass bit_thresh and inc == '!')
    seq2hmm[seq] = {model: [sstart, ssend, length, strand, score]} 
2) determine which model (archaea, bacteria, eukarya) the sequence most closely matches
    seq2hmm[seq] = [model, sstart, send, length, strand, score], [model2, sstart2, send2, length2, strand2, score2], ...]
3) identify regions that match to 23S (for best model)
4) mask internal regions that do not align to model
5) length threshold applies to aligned regions of 23S sequence
5) export 23S sequnece based on complete gene (including masked insertions)