def find_16S(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 16S (for best model)
4) mask internal regions that do not align to model
5) length threshold applies to aligned regions of 16S sequence
5) export 16S sequnece based on complete gene (including masked insertions)