Remove cluster with less than 10 sequences and
    loci with size smaller than 60%
435871
def _clean_cluster(list_c):
    
    global REMOVED
    init = len(list_c)
    list_c = {k: v for k, v in list_c.iteritems() if len(_get_seqs(v)) > parameters.min_seqs}
    logger.debug("_clean_cluster: number of clusters %s " % len(list_c.keys()))
    list_c = {k: _select_loci(v) for k, v in list_c.iteritems()}
    end = len(list_c)
    REMOVED += init - end
    return list_c
1060718
def _find_families_deprecated(clus_obj, min_seqs):
    
    logger.info("Creating meta-clusters based on shared sequences.")
    seen = defaultdict()
    metacluster = defaultdict(list)
    c_index = clus_obj.keys()
    meta_idx = 0
    with ProgressBar(maxval=len(c_index), redirect_stdout=True) as p:
        for itern, c in enumerate(c_index):
            p.update(itern)
            clus = clus_obj[c]
            if len(clus.idmembers.keys()) < min_seqs:
                del clus_obj[c]
                continue
            logger.debug("reading cluster %s" % c)
            logger.debug("loci2seq  %s" % clus.loci2seq)
            already_in, not_in = _get_seqs_from_cluster(clus.idmembers.keys(), seen)
            logger.debug("seen %s news %s" % (already_in, not_in))
            meta_idx += 1
            metacluster[meta_idx].append(c)
            seen.update(dict(zip(not_in, [meta_idx] * len(not_in))))
            if len(already_in) > 0:
                logger.debug("seen in %s" % already_in)
                for eindex in already_in:
                    for cluster in metacluster[eindex]:
                        metacluster[meta_idx].append(cluster)
                        prev_clus = clus_obj[cluster]
                        logger.debug("_find_families: prev %s current %s" % (eindex, clus.id))
                        # add current seqs to seen cluster
                        seqs_in = prev_clus.idmembers.keys()
                        seen.update(dict(zip(seqs_in, [meta_idx] * len(seqs_in))))
                        # for s_in_clus in prev_clus.idmembers:
                        #    seen[s_in_clus] = meta_idx
                    #    clus.idmembers[s_in_clus] = 1
                    # add current locus to seen cluster
                    # for loci in prev_clus.loci2seq:
                    #    logger.debug("adding %s" % loci)
                        # if not loci_old in current_clus[eindex].loci2seq:
                    #    clus.add_id_member(list(prev_clus.loci2seq[loci]), loci)
                    # logger.debug("loci %s" % clus.loci2seq.keys())
                    del metacluster[eindex]
                # clus_obj[c] = clus

                # logger.debug("num cluster %s" % len(clus_obj.keys()))
    logger.info("%s clusters merged" % len(metacluster))

    return metacluster, seen
1352231
private void removePairs(GvmCluster<S,K> cluster) {
		GvmClusterPair<S,K>[] pairs = cluster.pairs;
		for (int i = 0; i < bound-1; i++) {
			GvmClusterPair<S,K> pair = pairs[i];
			if (pair.c1.removed || pair.c2.removed) continue;
			this.pairs.remove(pair);
		}
		
	}
907931
def _remove_otiose(lst):
    """"""
    listtype = type([])
    while type(lst) == listtype and len(lst) == 1:
        lst = lst[0]

    return lst
545640
def _get_larger_chroms(ref_file):
    
    from scipy.cluster.vq import kmeans, vq
    all_sizes = []
    for c in ref.file_contigs(ref_file):
        all_sizes.append(float(c.size))
    all_sizes.sort()
    if len(all_sizes) > 5:
        # separate out smaller chromosomes and haplotypes with kmeans
        centroids, _ = kmeans(np.array(all_sizes), 2)
        idx, _ = vq(np.array(all_sizes), centroids)
        little_sizes = tz.first(tz.partitionby(lambda xs: xs[0], zip(idx, all_sizes)))
        little_sizes = [x[1] for x in little_sizes]
        # create one more cluster with the smaller, removing the haplotypes
        centroids2, _ = kmeans(np.array(little_sizes), 2)
        idx2, _ = vq(np.array(little_sizes), centroids2)
        little_sizes2 = tz.first(tz.partitionby(lambda xs: xs[0], zip(idx2, little_sizes)))
        little_sizes2 = [x[1] for x in little_sizes2]
        # get any chromosomes not in haplotype/random bin
        thresh = max(little_sizes2)
    else:
        thresh = 0
    larger_chroms = []
    for c in ref.file_contigs(ref_file):
        if c.size > thresh:
            larger_chroms.append(c.name)
    return larger_chroms
950573
def clean_inputs(data):
    
    if not utils.get_in(data, ("config", "algorithm", "variant_regions_orig")):
        data["config"]["algorithm"]["variant_regions_orig"] = dd.get_variant_regions(data)
    clean_vr = clean_file(dd.get_variant_regions(data), data, prefix="cleaned-")
    merged_vr = merge_overlaps(clean_vr, data)
    data["config"]["algorithm"]["variant_regions"] = clean_vr
    data["config"]["algorithm"]["variant_regions_merged"] = merged_vr

    if dd.get_coverage(data):
        if not utils.get_in(data, ("config", "algorithm", "coverage_orig")):
            data["config"]["algorithm"]["coverage_orig"] = dd.get_coverage(data)
        clean_cov_bed = clean_file(dd.get_coverage(data), data, prefix="cov-", simple=True)
        merged_cov_bed = merge_overlaps(clean_cov_bed, data)
        data["config"]["algorithm"]["coverage"] = clean_cov_bed
        data["config"]["algorithm"]["coverage_merged"] = merged_cov_bed

    if 'seq2c' in get_svcallers(data):
        seq2c_ready_bed = prep_seq2c_bed(data)
        if not seq2c_ready_bed:
            logger.warning("Can't run Seq2C without a svregions or variant_regions BED file")
        else:
            data["config"]["algorithm"]["seq2c_bed_ready"] = seq2c_ready_bed
    elif regions.get_sv_bed(data):
        dd.set_sv_regions(data, clean_file(regions.get_sv_bed(data), data, prefix="svregions-"))
    return data
112257
def fit_cluster(data, assignments, k, P_init, R_init, means):
    
    for c in range(k):
        if data[:,assignments==c].shape[1] == 0:
            _, assignments = kmeans_pp(data, k)
    genes, cells = data.shape
    nb_gene_indices = np.array([True for i in range(genes)])
    for c in range(k):
        c_data = data[:,assignments==c]
        nb_gene_indices = nb_gene_indices & find_nb_genes(c_data)
    for c in range(k):
        c_data = data[:,assignments==c]
        nb_genes = c_data[nb_gene_indices,:]
        poisson_genes = c_data[~nb_gene_indices, :]
        P_init[nb_gene_indices, c], R_init[nb_gene_indices, c] = nb_fit(nb_genes)
        means[~nb_gene_indices, c] = poisson_genes.mean(1)
    return nb_gene_indices
1323543
def run_all2(protgroup, memornot, subsequences, base_outdir,
            protgroup_dict, protein_feathers_dir, date, errfile, impute_counts=True,
            cutoff_num_proteins=0, core_only_genes=None,
            length_filter_pid=.8, remove_correlated_feats=True,
            force_rerun_counts=False, force_rerun_percentages=False, force_rerun_pca=False):
    """"""
    import ssbio.utils

    # Need to set multiprocessing limit for scipy/numpy stuff if parallelizing anything
    import os
    os.environ['OMP_NUM_THREADS'] = '1'

    # First, filter down the protein group to the membrane/nonmembrane definition
    prots_filtered_feathers = get_protein_feather_paths(protgroup=protgroup, memornot=memornot,
                                                        protgroup_dict=protgroup_dict,
                                                        protein_feathers_dir=protein_feathers_dir,
                                                        core_only_genes=core_only_genes)
    num_proteins = len(prots_filtered_feathers)
    if num_proteins <= cutoff_num_proteins:
        return

    # Make output directories
    protscale = 'proteome_unscaled'
    outdir_d0 = ssbio.utils.make_dir(op.join(base_outdir, protscale))
    outdir_d1 = ssbio.utils.make_dir(op.join(outdir_d0, '-'.join(memornot)))
    outdir_final = ssbio.utils.make_dir(op.join(outdir_d1, '-'.join(protgroup)))

    if impute_counts:
        big_strain_counts_df = get_proteome_counts_impute_missing(prots_filtered_feathers=prots_filtered_feathers,
                                                                  outpath=op.join(outdir_final,
                                                                                  '{}-subsequence_proteome_IMP.fthr'.format(
                                                                                      date)),
                                                                  length_filter_pid=length_filter_pid,
                                                                  force_rerun=force_rerun_counts)

        big_strain_percents_df = get_proteome_percentages(counts_df=big_strain_counts_df,
                                                          outpath=op.join(outdir_final,
                                                                          '{}-subsequence_proteome_perc_IMP.fthr'.format(
                                                                                  date)),
                                                          force_rerun=force_rerun_percentages)
        pca_pickle = op.join(outdir_final, '{}-subsequence_pca.pckl'.format(date))
    # Divide by totals to get percentages in a new dataframe
    else:
        try:
            big_strain_percents_df = get_proteome_correct_percentages(prots_filtered_feathers=prots_filtered_feathers,
                                                                      outpath=op.join(outdir_final,
                                                                                      '{}-subsequence_proteome_perc_AVG.fthr'.format(
                                                                                              date)),
                                                                      length_filter_pid=length_filter_pid,
                                                                      force_rerun=force_rerun_percentages)
            pca_pickle = op.join(outdir_final, '{}-subsequence_pca_AVG.pckl'.format(date))
        except:
            with open(errfile, "a") as myfile:
                myfile.write('PERCENTAGES ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return



    if ssbio.utils.force_rerun(flag=force_rerun_pca, outfile=pca_pickle):

        # Then, get filters for rows of the loaded feathers for interested subsequences
        keep_subsequences = get_interested_subsequences(subsequences=subsequences)

        # Some numbers: number of features
        num_feats = len(big_strain_percents_df)

        # Make an unwieldy title
        big_title = 'LOC={0}; PROTGROUP={1};\n' \
                    'NUMPROTS={2}; NUMFEATS={3}'.format('-'.join(memornot),
                                                                        '-'.join(protgroup),
                                                                        num_proteins,
                                                                        num_feats)

        # Run PCA and make plots
        runner = PCAMultiROS(features_df=big_strain_percents_df, observations_df=pd.DataFrame(), plot_title=big_title)
        try:
            runner.clean_data(keep_features=keep_subsequences, remove_correlated_feats=remove_correlated_feats)
        except:
            with open(errfile, "a") as myfile:
                myfile.write(
                    'CLEAN ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return
        # try:
        runner.run_pca()
        # except:
        #     with open(errfile, "a") as myfile:
        #         myfile.write(
        #             'PCA ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
        #     return
        with open(pca_pickle, 'wb') as f:
            pickle.dump(runner, f)
    else:
        with open(pca_pickle, 'rb') as f:
            runner = pickle.load(f)
1101118
def _select_loci(c):
    """"""
    loci_len = {k: len(v) for k, v in c.loci2seq.iteritems()}
    logger.debug("_select_loci: number of loci %s" % len(c.loci2seq.keys()))
    loci_len_sort = sorted(loci_len.iteritems(), key=operator.itemgetter(1), reverse=True)
    max_size = loci_len_sort[0][1]
    logger.debug("_select_loci: max size %s" % max_size)
    loci_clean = {locus: c.loci2seq[locus] for locus, size in loci_len_sort if size > 0.8 * max_size}
    c.loci2seq = loci_clean
    removed = list(set(c.idmembers.keys()) - set(_get_seqs(c)))
    c.add_id_member(removed, loci_len_sort[0][0])
    logger.debug("_select_loci: number of loci %s after cleaning" % len(c.loci2seq.keys()))
    return c
62850
def cluster_seqs(seqs,
                 neighbor_join=False,
                 params={},
                 add_seq_names=True,
                 WorkingDir=tempfile.gettempdir(),
                 SuppressStderr=None,
                 SuppressStdout=None,
                 max_chars=1000000,
                 max_hours=1.0,
                 constructor=PhyloNode,
                 clean_up=True
                 ):
    
    num_seqs = len(seqs)
    if num_seqs < 2:
        raise ValueError, "Muscle requres 2 or more sequences to cluster."


    num_chars = sum(map(len, seqs))
    if num_chars > max_chars:
        params["-maxiters"] = 2
        params["-diags1"] = True
        params["-sv"] = True
        #params["-distance1"] = "kmer6_6"
        #params["-distance1"] = "kmer20_3"
        #params["-distance1"] = "kbit20_3"
        print "lots of chars, using fast align", num_chars


    params["-maxhours"] = max_hours
    #params["-maxiters"] = 10

    #cluster_type = "upgmb"
    #if neighbor_join:
    #    cluster_type = "neighborjoining"

    params["-clusteronly"] = True
    params["-tree1"] = get_tmp_filename(WorkingDir)

    muscle_res = muscle_seqs(seqs,
                 params=params,
                 add_seq_names=add_seq_names,
                 WorkingDir=WorkingDir,
                 SuppressStderr=SuppressStderr,
                 SuppressStdout=SuppressStdout)

    tree = DndParser(muscle_res["Tree1Out"], constructor=constructor)

    if clean_up:
        muscle_res.cleanUp()
    return tree
114473
def _remove_end_gaps(sequence):
    '''

    '''
    # Count terminal blank sequences
    def count_end_gaps(seq):
        gap = coral.DNA('-')
        count = 0
        for base in seq:
            if base == gap:
                count += 1
            else:
                break

        return count

    top_left = count_end_gaps(sequence.top)
    top_right = count_end_gaps(reversed(sequence.top))
    bottom_left = count_end_gaps(reversed(sequence.bottom))
    bottom_right = count_end_gaps(sequence.bottom)

    # Trim sequence
    left_index = min(top_left, bottom_left)
    right_index = len(sequence) - min(top_right, bottom_right)

    return sequence[left_index:right_index]
983725
def fake_run(self):
        ''''''
        clusters = {}
        used_names = set()
        seq_reader = pyfastaq.sequences.file_reader(self.infile)

        for seq in seq_reader:
            if seq.id in used_names:
                raise Error('Sequence name "' + seq.id + '" not unique. Cannot continue')

            clusters[str(len(clusters) + self.min_cluster_number)] = {seq.id}
            used_names.add(seq.id)

        return clusters
599145
def get_sequences_from_cluster(c1, c2, data):
    """"""
    seqs1 = data[c1]['seqs']
    seqs2 = data[c2]['seqs']
    seqs = list(set(seqs1 + seqs2))
    names = []
    for s in seqs:
        if s in seqs1 and s in seqs2:
            names.append("both")
        elif s in seqs1:
            names.append(c1)
        else:
            names.append(c2)
    return seqs, names
677102
def scan(self, seqs, nreport=100, scan_rc=True, normalize=False):
        

        if not self.threshold:
            sys.stderr.write(
                "Using default threshold of 0.95. "
                "This is likely not optimal!\n"
                )
            self.set_threshold(threshold=0.95)

        seqs = as_fasta(seqs, genome=self.genome)
           
        it = self._scan_sequences(seqs.seqs, 
                    nreport, scan_rc)
       
        if normalize:
            if len(self.meanstd) == 0:
                self.set_meanstd()
            mean_std = [self.meanstd.get(m_id) for m_id in self.motif_ids]
            means = [x[0] for x in  mean_std]
            stds = [x[1] for x in  mean_std]


        for result in it:
            if normalize:
                zresult = [] 
                for i,mrow in enumerate(result):
                    mrow = [((x[0] - means[i]) / stds[i], x[1], x[2]) for x in mrow]
                    zresult.append(mrow)
                yield zresult
            else:
                yield result
225571
def _truncate_seq_pair(self, tokens_a, tokens_b, max_length):
        """"""
        # This is a simple heuristic which will always truncate the longer sequence
        # one token at a time. This makes more sense than truncating an equal percent
        # of tokens from each, since if one sequence is very short then each token
        # that's truncated likely contains more information than a longer sequence.
        while True:
            total_length = len(tokens_a) + len(tokens_b)
            if total_length <= max_length:
                break
            if len(tokens_a) > len(tokens_b):
                tokens_a.pop()
            else:
                tokens_b.pop()
592952
def _regroup(self):
        
        # clean up old maps
        for j in range(self.nout):
            self.inv_map[j] = []

        # find smallest divergence between input component i
        # and output component j of the cluster mixture density
        for i in range(self.nin):
            self.min_kl[i] = np.inf
            j_min = None
            for j in range(self.nout):
                kl = kullback_leibler(self.f.components[i], self.g.components[j])
                if kl < self.min_kl[i]:
                    self.min_kl[i] = kl
                    j_min = j
            assert j_min is not None
            self.inv_map[j_min].append(i)
753669
def remove_variants(self, variants):
        """"""
        chroms = set([i.chrom for i in variants])
        for chrom in chroms:
            if self.append_chromosome:
                chrom = 'chr%s' % chrom
            to_delete = [pos for pos in self.positions[chrom] if pos in variants]
            for pos in to_delete:
                del self.positions[chrom][pos]