Prepare alignment for cluster detection.
179517
def _clean_alignment(args):
    
    logger.info("Clean bam file with highly repetitive reads with low counts. sum(counts)/n_hits > 1%")
    bam_file, seq_obj = clean_bam_file(args.afile, args.mask)
    logger.info("Using %s file" % bam_file)
    detect_complexity(bam_file, args.ref, args.out)
    return bam_file, seq_obj
262661
static void pruneIncorrectSize(List<List<EllipsesIntoClusters.Node>> clusters, int N) {
		// prune clusters which can't be a member calibration target
		for (int i = clusters.size()-1; i >= 0; i--) {
			if( clusters.get(i).size() != N ) {
				clusters.remove(i);
			}
		}
	}
227310
def estimate_map(interface, state, label, inp):
    """"""
    out = interface.output(0)
    centers = {}
    for row in inp:
        row = row.strip().split(state["delimiter"])
        if len(row) > 1:
            x = [(0 if row[i] in state["missing_vals"] else float(row[i])) for i in state["X_indices"]]
            cluster = min((state['dist'](c, x), i) for i, c in state['centers'])[1]
            vertex = state['create'](x, 1.0)
            centers[cluster] = vertex if cluster not in centers else state["update"](centers[cluster], vertex)

    for cluster, values in centers.iteritems():
        out.add(cluster, values)
1196532
def init_pp(X, n_clusters, random_state):
    
    x_squared_norms = row_norms(X, squared=True).compute()
    logger.info("Initializing with k-means++")
    with _timer("initialization of %2d centers" % n_clusters, _logger=logger):
        centers = sk_k_means._k_init(
            X, n_clusters, random_state=random_state, x_squared_norms=x_squared_norms
        )

    return centers
1241492
def eeg_microstates_clustering(data, n_microstates=4, clustering_method="kmeans", n_jobs=1, n_init=25, occurence_rejection_treshold=0.05, max_refitting=5, verbose=True):
    
    # Create training set
    training_set = data.copy()

    if verbose is True:
        print("- Initializing the clustering algorithm...")
    if clustering_method == "kmeans":
        algorithm = sklearn.cluster.KMeans(init='k-means++', n_clusters=n_microstates, n_init=n_init, n_jobs=n_jobs)
    elif clustering_method == "spectral":
        algorithm = sklearn.cluster.SpectralClustering(n_clusters=n_microstates, n_init=n_init, n_jobs=n_jobs)
    elif clustering_method == "agglom":
        algorithm = sklearn.cluster.AgglomerativeClustering(n_clusters=n_microstates, linkage="complete")
    elif clustering_method == "dbscan":
        algorithm = sklearn.cluster.DBSCAN(min_samples=100)
    elif clustering_method == "affinity":
        algorithm = sklearn.cluster.AffinityPropagation(damping=0.5)
    else:
        print("NeuroKit Error: eeg_microstates(): clustering_method must be 'kmeans', 'spectral', 'dbscan', 'affinity' or 'agglom'")


    refitting = 0  # Initialize the number of refittings
    good_fit_achieved = False
    while good_fit_achieved is False:
        good_fit_achieved = True
        if verbose is True:
            print("- Fitting the classifier...")
        # Fit the algorithm
        algorithm.fit(training_set)

        if verbose is True:
            print("- Clustering back the initial data...")
        # Predict the more likely cluster for each observation
        predicted = algorithm.fit_predict(training_set)

        if verbose is True:
            print("- Check for abnormalities...")
        # Check for abnormalities and prune the training set until none found
        occurences = dict(collections.Counter(predicted))
        masks = [np.array([True]*len(training_set))]
        for microstate in occurences:
            # is the frequency of one microstate inferior to a treshold
            if occurences[microstate] < len(data)*occurence_rejection_treshold:
                good_fit_achieved = False
                refitting += 1  # Increment the refitting
                print("NeuroKit Warning: eeg_microstates(): detected some outliers: refitting the classifier (n=" + str(refitting) + ").")
                masks.append(predicted!=microstate)
        mask = np.all(masks, axis=0)
        training_set = training_set[mask]

    return(algorithm)
50007
def initialize_means_weights(data, clusters, init_means=None, init_weights=None, initialization='tsvd', max_assign_weight=0.75):
    
    genes, cells = data.shape
    if init_means is None:
        if init_weights is not None:
            if len(init_weights.shape)==1:
                means = initialize_means(data, init_weights, clusters)
            else:
                means = initialize_means(data, init_weights.argmax(0),
                        clusters, max_assign_weight=max_assign_weight)
        elif initialization=='cluster':
            assignments, means = poisson_cluster(data, clusters)
            if init_weights is None:
                init_weights = initialize_from_assignments(assignments, clusters,
                        max_assign_weight=max_assign_weight)
        elif initialization=='kmpp':
            means, assignments = kmeans_pp(data, clusters)
        elif initialization=='km':
            km = KMeans(clusters)
            assignments = km.fit_predict(log1p(cell_normalize(data)).T)
            init_weights = initialize_from_assignments(assignments, clusters,
                    max_assign_weight)
            means = initialize_means(data, assignments, clusters)
        elif initialization=='tsvd':
            n_components = min(50, genes-1)
            #tsvd = TruncatedSVD(min(50, genes-1))
            km = KMeans(clusters)
            # remove dependence on sklearn tsvd b/c it has a bug that
            # prevents it from working properly on long inputs 
            # if num elements > 2**31
            #data_reduced = tsvd.fit_transform(log1p(cell_normalize(data)).T)
            U, Sigma, VT = randomized_svd(log1p(cell_normalize(data)).T,
                    n_components)
            data_reduced = U*Sigma
            assignments = km.fit_predict(data_reduced)
            init_weights = initialize_from_assignments(assignments, clusters,
                    max_assign_weight)
            means = initialize_means(data, assignments, clusters)
        elif initialization == 'random' or initialization == 'rand':
            # choose k random cells and set means to those
            selected_cells = np.random.choice(range(cells), size=clusters,
                    replace=False)
            means = data[:, selected_cells]
            if sparse.issparse(means):
                means = means.toarray()
    else:
        means = init_means.copy()
    means = means.astype(float)
    if init_weights is None:
        if init_means is not None:
            if initialization == 'cluster':
                assignments, means = poisson_cluster(data, clusters,
                        init=init_means, max_iters=1)
                w_init = initialize_from_assignments(assignments, clusters,
                        max_assign_weight)
            elif initialization == 'km':
                km = KMeans(clusters, init=log1p(init_means.T), max_iter=1)
                assignments = km.fit_predict(log1p(cell_normalize(data)).T)
                w_init = initialize_from_assignments(assignments, clusters,
                        max_assign_weight)
            else:
                w_init = np.random.random((clusters, cells))
                w_init = w_init/w_init.sum(0)
        else:
            w_init = np.random.random((clusters, cells))
            w_init = w_init/w_init.sum(0)
    else:
        if len(init_weights.shape)==1:
            init_weights = initialize_from_assignments(init_weights, clusters,
                    max_assign_weight)
        w_init = init_weights.copy()
    return means, w_init
1038266
def align(fastq_file, pair_file, index_dir, names, align_dir, data):
    
    umi_ext = "-cumi" if "umi_bam" in data else ""
    out_file = os.path.join(align_dir, "{0}-sort{1}.bam".format(dd.get_sample_name(data), umi_ext))
    num_cores = data["config"]["algorithm"].get("num_cores", 1)
    rg_info = "rgid={rg} rgpl={pl} rgpu={pu} rgsm={sample}".format(**names)
    pair_file = pair_file if pair_file else ""
    final_file = None
    if data.get("align_split"):
        # BBMap does not accept input fastq streams
        raise ValueError("bbmap is not compatible with alignment splitting, set `align_split: false`")
    pair_arg = "in2=%s" % pair_file if pair_file else ""
    if not utils.file_exists(out_file) and (final_file is None or not utils.file_exists(final_file)):
        with postalign.tobam_cl(data, out_file, pair_file != "") as (tobam_cl, tx_out_file):
            if index_dir.endswith(("/ref", "/ref/")):
                index_dir = os.path.dirname(index_dir)
            # sam=1.3 required for compatibility with strelka2
            cmd = ("bbmap.sh sam=1.3 mdtag=t {rg_info} path={index_dir} in1={fastq_file} "
                   "{pair_arg} out=stdout.sam | ")
            do.run(cmd.format(**locals()) + tobam_cl, "bbmap alignment: %s" % dd.get_sample_name(data))
    data["work_bam"] = out_file
    return data
1165550
def get_cluster_assignment(self):
        """"""
        plan = self.get_cluster_plan()
        assignment = {}
        for elem in plan['partitions']:
            assignment[
                (elem['topic'], elem['partition'])
            ] = elem['replicas']

        return assignment
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
964052
def read_cluster(data, id=1):
    """"""
    cl = cluster(1)

    # seqs = [s.values()[0] for s in data['seqs']]
    names = [s.keys()[0] for s in data['seqs']]
    cl.add_id_member(names, 1)
    freq = defaultdict()
    [freq.update({s.keys()[0]: s.values()[0]}) for s in data['freq']]
659623
public static int nvgraphAnalyzeClustering(
        nvgraphHandle handle, 
        nvgraphGraphDescr graph_descr, 
        long weight_index, 
        int n_clusters, 
        Pointer clustering, 
        int metric, 
        Pointer score)
    {
        return checkResult(nvgraphAnalyzeClusteringNative(handle, graph_descr, weight_index, n_clusters, clustering, metric, score));
    }
1202010
protected void computeClusterDistance() {
		for( int i = 0; i < pixels.size; i++ ) {
			pixels.data[i].reset();
		}

		for( int i = 0; i < clusters.size && !stopRequested; i++ ) {
			Cluster c = clusters.data[i];

			// compute search bounds
			int centerX = (int)(c.x + 0.5f);
			int centerY = (int)(c.y + 0.5f);

			int x0 = centerX - gridInterval; int x1 = centerX + gridInterval + 1;
			int y0 = centerY - gridInterval; int y1 = centerY + gridInterval + 1;

			if( x0 < 0 ) x0 = 0;
			if( y0 < 0 ) y0 = 0;
			if( x1 > input.width ) x1 = input.width;
			if( y1 > input.height ) y1 = input.height;

			for( int y = y0; y < y1; y++ ) {
				int indexPixel = y*input.width + x0;
				int indexInput = input.startIndex + y*input.stride + x0;

				int dy = y-centerY;

				for( int x = x0; x < x1; x++ ) {
					int dx = x-centerX;

					float distanceColor = colorDistance(c.color,indexInput++);
					float distanceSpacial = dx*dx + dy*dy;
					pixels.data[indexPixel++].add(c,distanceColor + adjustSpacial*distanceSpacial);
				}
			}
		}
	}
176821
protected void align() {
		if (!isReady()) {
			return;
		}

		long timeStart = System.nanoTime();

		int[] dim = getScoreMatrixDimensions();
		if (storingScoreMatrix) {
			scores = new int[dim[0]][dim[1]][dim[2]];
		} else {
			scores = new int[dim[0]][][];
			scores[0] = new int[dim[1]][dim[2]];
			scores[1] = new int[dim[1]][dim[2]];
		}
		boolean linear = (gapPenalty.getType() == GapPenalty.Type.LINEAR);
		Last[][][] traceback = new Last[dim[0]][][];
		List<Step> sx = new ArrayList<Step>(), sy = new ArrayList<Step>();

		if (!local) {
			xyMax = new int[] { dim[0] - 1, dim[1] - 1 };
			xyStart = new int[] { 0, 0 };
			score = 0;
			List<Subproblem> problems = Subproblem.getSubproblems(anchors, xyMax[0], xyMax[1]);
			assert problems.size() == anchors.size() + 1;
			for (int i = 0; i < problems.size(); i++) {
				Subproblem subproblem = problems.get(i);
				for (int x = subproblem.getQueryStartIndex(); x <= subproblem.getQueryEndIndex(); x++) {

					traceback[x] =

					  linear ?
						setScoreVector(x, subproblem, gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x, subproblem), storingScoreMatrix, scores) :

						setScoreVector(x, subproblem, gapPenalty.getOpenPenalty(), gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x, subproblem), storingScoreMatrix, scores);
				}
			}
			setSteps(traceback, scores, sx, sy);
			score = Integer.MIN_VALUE;
			int[] finalScore = scores[xyMax[0]][xyMax[1]];
			for (int z = 0; z < finalScore.length; z++) {
				score = Math.max(score, finalScore[z]);
			}
		} else {
			for (int x = 0; x < dim[0]; x++) {

				traceback[x] =

				  linear ?
					setScoreVector(x, gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x), storingScoreMatrix, scores, xyMax, score) :

					setScoreVector(x, gapPenalty.getOpenPenalty(), gapPenalty.getExtensionPenalty(), getSubstitutionScoreVector(x), storingScoreMatrix, scores, xyMax, score);

				if (xyMax[0] == x) {
					score = scores[x][xyMax[1]][0];
				}
			}
			xyStart = local ? setSteps(traceback, xyMax, sx, sy) : setSteps(traceback, scores, sx, sy);
		}

		setProfile(sx, sy);

		if (!storingScoreMatrix) {
			scores = null;
		}

		time = System.nanoTime() - timeStart;
	}
687463
def _split_cluster(c, pairs, n):
    """"""
    old = c[p[0]]
    new = c[p[1]]
    new_c = cluster(n)
    common = set(_get_seqs(old)).intersection(_get_seqs(new))
    for idl in old.loci2seq:
        in_common = list(set(common).intersection(old.loci2seq[idl]))
        if len(in_common) > 0:
            logger.debug("_split_cluster: in_common %s with pair 1" % (len(in_common)))
            new_c.add_id_member(in_common, idl)
            old.loci2seq[idl] = list(set(old.loci2seq[idl]) - set(common))
            logger.debug("_split_cluster: len old %s with pair 1" % (len(old.loci2seq)))
    for idl in new.loci2seq:
        in_common = list(set(common).intersection(new.loci2seq[idl]))
        if len(in_common) > 0:
            logger.debug("_split_cluster: in_common %s with pair 2" % (len(in_common)))
            new_c.add_id_member(in_common, idl)
            new.loci2seq[idl] = list(set(new.loci2seq[idl]) - set(common))
            logger.debug("_split_cluster: len old %s with pair 2" % (len(new.loci2seq)))
    old.update()
    new.update()
    old.loci2seq = {k: v for k, v in old.loci2seq.iteritems() if len(v) > 0}
    new.loci2seq = {k: v for k, v in new.loci2seq.iteritems() if len(v) > 0}
    c[n] = new
    c[p[0]] = old
    c[p[1]] = new
    return c
12299
def map_alignment(self, prot_seq, nucl_seq):
        
        if prot_seq.id != nucl_seq.id:
            logging.warning(
                'ID mismatch: %s != %s. Are the sequences in the same order?',
                prot_seq.id, nucl_seq.id)

        # Ungap nucleotides
        codons = batch(str(nucl_seq.seq.ungap('-')), 3)
        codons = [''.join(i) for i in codons]
        codon_iter = iter(codons)

        ungapped_prot = str(prot_seq.seq).replace('-', '')

        if len(ungapped_prot) != len(codons):
            table = self.translation_table.forward_table
            prot_str = ' '.join(' ' + p + ' ' for p in ungapped_prot)
            codon_str = ' '.join(codons)
            trans_str = ' '.join(' ' + table.get(codon, 'X') + ' '
                                 for codon in codons)
            raise ValueError("""Length of codon sequence ({0}) does not match \
length of protein sequence ({1}) for {2}
Protein:       {3}
Codons:        {4}
Trans. Codons: {5}""".format(len(codons), len(ungapped_prot), nucl_seq.id, prot_str,
                            codon_str, trans_str))

        try:
            nucl_align = ['---' if p == '-' else next(codon_iter)
                          for p in str(prot_seq.seq)]
        except StopIteration:
            assert False  # Should be checked above

        result = SeqRecord(Seq(''.join(nucl_align)), id=nucl_seq.id,
                           description=nucl_seq.description)

        # Validate
        self._validate_translation(prot_seq.seq.upper(), result.seq.upper())

        return result
714849
def fiGroupFromLinkageGroup(matchArr, arrPos, groupId, timeKey, massKey):
    
    fgi = Fgi(groupId)
    matchArr['isAnnotated'][arrPos]

    minT = numpy.min(matchArr[timeKey][arrPos])
    maxT = numpy.max(matchArr[timeKey][arrPos])
    minM = numpy.min(matchArr[massKey][arrPos])
    maxM = numpy.max(matchArr[massKey][arrPos])

    fgi.isValid = True
    fgi.specfiles = matchArr['specfile'][arrPos].tolist()
    fgi.featureIds = matchArr['id'][arrPos].tolist()
    fgi.isAnnotated = numpy.any(matchArr['isAnnotated'][arrPos])
    fgi.coordinates = ((minT, maxT), (minM, maxM))
    #fgi.clusterType = clusterType

    return fgi
883936
def cluster_labels_B(hdf5_file, s_reduced, lock, I, ii, iix, rows_slice):
    

    with Worker.hdf5_lock:
        with tables.open_file(hdf5_file, 'r+') as fileh:
            S = fileh.root.aff_prop_group.similarities
            s = S[rows_slice, ...]
           
    s = s[:, ii]
    s = s[iix[rows_slice]]
    
    with lock:                
        s_reduced += s[:].sum(axis = 0)
        
    del s