<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def protein_permutation(graph_score, num_codons_obs, context_counts, context_to_mut, seq_context, gene_seq, gene_graph, num_permutations=10000, stop_criteria=100, pseudo_count=0):
 """Performs null-simulations for position-based mutation statistics in a single gene. Parameters graph_score : float clustering score for observed data num_codons_obs : int number of codons with missense mutation in observed data context_counts : pd.Series number of mutations for each context context_to_mut : dict dictionary mapping nucleotide context to a list of observed somatic base changes. seq_context : SequenceContext Sequence context for the entire gene sequence (regardless of where mutations occur). The nucleotide contexts are identified at positions along the gene. gene_seq : GeneSequence Sequence of gene of interest num_permutations : int, default: 10000 number of permutations to create for null stop_criteria : int stop after stop_criteria iterations are more significant then the observed statistic. Returns ------- protein_pval : float p-value for clustering in neighbor graph constructure from protein structures """

    # get contexts and somatic base
    mycontexts = context_counts.index.tolist()
    somatic_base = [base
                    for one_context in mycontexts
                    for base in context_to_mut[one_context]]

    # get random positions determined by sequence context
    tmp_contxt_pos = seq_context.random_pos(context_counts.iteritems(),
                                            num_permutations)
    tmp_mut_pos = np.hstack(pos_array for base, pos_array in tmp_contxt_pos)

    # calculate position-based statistics as a result of random positions
    null_graph_entropy_ct = 0
    coverage_list = []
    num_mut_list = []
    graph_entropy_list = []
    for i, row in enumerate(tmp_mut_pos):
        # calculate the expected value of the relative increase in coverage
        if i == stop_criteria-1:
            rel_inc = [coverage_list[k] / float(num_mut_list[k])
                       for k in range(stop_criteria-1)
                       if coverage_list[k]]
            exp_rel_inc = np.mean(rel_inc)

            # calculate observed statistic
            if num_codons_obs:
                obs_stat = graph_score / np.log2(exp_rel_inc*num_codons_obs)
            else:
                obs_stat = 1.0

            # calculate statistics for simulated data
            sim_stat_list = [ent / np.log2(exp_rel_inc*num_mut_list[l])
                             for l, ent in enumerate(graph_entropy_list)]
            null_graph_entropy_ct = len([s for s in sim_stat_list
                                         if s-utils.epsilon <= obs_stat])

        # get info about mutations
        tmp_mut_info = mc.get_aa_mut_info(row,
                                          somatic_base,
                                          gene_seq)

        # calculate position info
        tmp_tuple = cutils.calc_pos_info(tmp_mut_info['Codon Pos'],
                                         tmp_mut_info['Reference AA'],
                                         tmp_mut_info['Somatic AA'],
                                         pseudo_count=pseudo_count,
                                         is_obs=0)
        _, _, _, tmp_pos_ct = tmp_tuple

        # record num of mut codons
        if i < stop_criteria-1:
            tmp_num_mut_codons = len(tmp_pos_ct)
            num_mut_list.append(tmp_num_mut_codons)

        # get entropy on graph-smoothed probability distribution
        tmp_graph_entropy, tmp_coverage = scores.compute_ng_stat(gene_graph, tmp_pos_ct)

        # record the "coverage" in the graph
        if i < stop_criteria-1:
            coverage_list.append(tmp_coverage)
            graph_entropy_list.append(tmp_graph_entropy)

        # update empirical null distribution counts
        if i >= stop_criteria:
            #if tmp_graph_entropy-utils.epsilon <= graph_score:
            if tmp_num_mut_codons:
                sim_stat = tmp_graph_entropy / np.log2(exp_rel_inc*tmp_num_mut_codons)
            else:
                sim_stat = 1.0

            # add count
            if sim_stat-utils.epsilon <= obs_stat:
                null_graph_entropy_ct += 1

        # stop iterations if reached sufficient precision
        if null_graph_entropy_ct >= stop_criteria:
            break

    # calculate p-value from empirical null-distribution
    protein_pval = float(null_graph_entropy_ct) / (i+1)

    return protein_pval, obs_stat