def hotmaps_permutation(obs_stat,
                        context_counts,
                        context_to_mut,
                        seq_context,
                        gene_seq,
                        window,
                        num_permutations=10000,
                        stop_criteria=100,
                        max_batch=25000,
                        null_save_path=None):
    """Performs null-permutations for position-based mutation statistics
    in a single gene.

    Parameters
    ----------
    obs_stat : dict
        dictionary mapping codons to the sum of mutations in a window
    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
    window : int
        Number of codons to the left/right of a mutated position to consider
        in the window
    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.
    max_batch : int
        maximum number of whole gene simulations to do at once.
        For large number of simulations holding a matrix of M x N,
        where M is the number of mutations and N is the number of simulations,
        can get quite large.
    null_save_path : str or None
        File path to save null distribution. If None, don't save it.

    Returns
    -------
    pvals : dict
        Maps mutated codon position to the calculated p-value
    """
    # 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]]

    # calculate the # of batches for simulations
    max_batch = min(num_permutations, max_batch)
    num_batches = num_permutations // max_batch
    remainder = num_permutations % max_batch
    batch_sizes = [max_batch] * num_batches
    if remainder:
        batch_sizes += [remainder]

    # figure out which position has highest value
    max_key = {w: max(obs_stat[w], key=(lambda key: obs_stat[w][key]))
               for w in window}

    # setup null dist counts
    null_cts = {w: {k: 0 for k in obs_stat[w]}
                for w in window }

    # empirical null distribution (saved if file path provided)
    empirical_null = {w: {} for w in window}

    num_sim = 0 # number of simulations
    for j, batch_size in enumerate(batch_sizes):
        # stop iterations if reached sufficient precision
        # stop iterations if reached sufficient precision
        stop_flag = [(null_cts[w][max_key[w]]>=stop_criteria)
                      for w in window]
        if all(stop_flag):
            break
        #if null_cts[max_key] >= stop_criteria:
            #break

        # get random positions determined by sequence context
        tmp_contxt_pos = seq_context.random_pos(context_counts.iteritems(),
                                                batch_size)
        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
        for i, row in enumerate(tmp_mut_pos):
            # get info about mutations
            tmp_mut_info = mc.get_aa_mut_info(row,
                                              somatic_base,
                                              gene_seq)

            # calculate position info
            tmp_pos, tmp_sim = utils.calc_windowed_sum(tmp_mut_info['Codon Pos'],
                                                 tmp_mut_info['Reference AA'],
                                                 tmp_mut_info['Somatic AA'],
                                                 window)

            # update the counts when the empirical null passes the observed
            for tmp_w in tmp_sim:
                for tmp_key in tmp_sim[tmp_w]:
                    # get mutation count for simulation
                    val = tmp_sim[tmp_w][tmp_key]

                    # add to empirical null distribution
                    empirical_null[tmp_w].setdefault(val, 0)
                    empirical_null[tmp_w][val] += 1

                    # update counts used for p-value
                    for key in null_cts[tmp_w]:
                        if val >= obs_stat[tmp_w][key]:
                            null_cts[tmp_w][key] += 1

            # update the number of simulations
            num_sim += len(tmp_pos)

            # stop iterations if reached sufficient precision
            stop_flag = [(null_cts[w][max_key[w]]>=stop_criteria)
                         for w in window]
            if all(stop_flag):
                break

    # calculate p-value from empirical null-distribution
    pvals = {w: {k: float(null_cts[w][k]) / (num_sim) for k in obs_stat[w]}
             for w in window}

    # save empirical distribution
    if null_save_path:
        for w in window:
            # create null distribution
            output = [['mutation_count', 'p-value']]
            sorted_cts = sorted(empirical_null[w].keys())
            tmp_sum = 0
            for i in range(len(sorted_cts)):
                tmp_sum += empirical_null[w][sorted_cts[-(i+1)]]
                tmp_pval = tmp_sum / float(num_sim)
                output.append([sorted_cts[-(i+1)], tmp_pval])
            # save output
            with open(null_save_path.format(w), 'w') as handle:
                mywriter = csv.writer(handle, delimiter='\t', lineterminator='\n')
                mywriter.writerows(output)

    return pvals
def hotmaps_permutation ( obs_stat , context_counts , context_to_mut , seq_context , gene_seq , window , num_permutations = 10000 , stop_criteria = 100 , max_batch = 25000 , null_save_path = None ) : mycontexts = context_counts . index . tolist ( ) somatic_base = [ base for one_context in mycontexts for base in context_to_mut [ one_context ] ] max_batch = min ( num_permutations , max_batch ) num_batches = num_permutations // max_batch remainder = num_permutations % max_batch batch_sizes = [ max_batch ] * num_batches if remainder : batch_sizes += [ remainder ] max_key = { w : max ( obs_stat [ w ] , key = ( lambda key : obs_stat [ w ] [ key ] ) ) for w in window } null_cts = { w : { k : 0 for k in obs_stat [ w ] } for w in window } empirical_null = { w : { } for w in window } num_sim = 0 for j , batch_size in enumerate ( batch_sizes ) : stop_flag = [ ( null_cts [ w ] [ max_key [ w ] ] >= stop_criteria ) for w in window ] if all ( stop_flag ) : break tmp_contxt_pos = seq_context . random_pos ( context_counts . iteritems ( ) , batch_size ) tmp_mut_pos = np . hstack ( pos_array for base , pos_array in tmp_contxt_pos ) for i , row in enumerate ( tmp_mut_pos ) : tmp_mut_info = mc . get_aa_mut_info ( row , somatic_base , gene_seq ) tmp_pos , tmp_sim = utils . calc_windowed_sum ( tmp_mut_info [ 'Codon Pos' ] , tmp_mut_info [ 'Reference AA' ] , tmp_mut_info [ 'Somatic AA' ] , window ) for tmp_w in tmp_sim : for tmp_key in tmp_sim [ tmp_w ] : val = tmp_sim [ tmp_w ] [ tmp_key ] empirical_null [ tmp_w ] . setdefault ( val , 0 ) empirical_null [ tmp_w ] [ val ] += 1 for key in null_cts [ tmp_w ] : if val >= obs_stat [ tmp_w ] [ key ] : null_cts [ tmp_w ] [ key ] += 1 num_sim += len ( tmp_pos ) stop_flag = [ ( null_cts [ w ] [ max_key [ w ] ] >= stop_criteria ) for w in window ] if all ( stop_flag ) : break pvals = { w : { k : float ( null_cts [ w ] [ k ] ) / ( num_sim ) for k in obs_stat [ w ] } for w in window } if null_save_path : for w in window : output = [ [ 'mutation_count' , 'p-value' ] ] sorted_cts = sorted ( empirical_null [ w ] . keys ( ) ) tmp_sum = 0 for i in range ( len ( sorted_cts ) ) : tmp_sum += empirical_null [ w ] [ sorted_cts [ - ( i + 1 ) ] ] tmp_pval = tmp_sum / float ( num_sim ) output . append ( [ sorted_cts [ - ( i + 1 ) ] , tmp_pval ] ) with open ( null_save_path . format ( w ) , 'w' ) as handle : mywriter = csv . writer ( handle , delimiter = '\t' , lineterminator = '\n' ) mywriter . writerows ( output ) return pvals
Performs null-permutations for position-based mutation statistics
    in a single gene.

    Parameters
    ----------
    obs_stat : dict
        dictionary mapping codons to the sum of mutations in a window
    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
    window : int
        Number of codons to the left/right of a mutated position to consider
        in the window
    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.
    max_batch : int
        maximum number of whole gene simulations to do at once.
        For large number of simulations holding a matrix of M x N,
        where M is the number of mutations and N is the number of simulations,
        can get quite large.
    null_save_path : str or None
        File path to save null distribution. If None, don't save it.

    Returns
    -------
    pvals : dict
        Maps mutated codon position to the calculated p-value