m.bestseqs(,thresh=None)  -- Return all k-mers that match motif with a score >= thresh
190009
def bestseqs(self,thresh=None):
        
        if not thresh:
            if self._bestseqs:
                return self._bestseqs
        if not thresh: thresh = 0.8 * self.maxscore
        self._bestseqs = bestseqs(self,thresh)
        return self._bestseqs
513323
def cluster(args):
    

    args = _check_args(args)
    read_stats_file = op.join(args.dir_out, "read_stats.tsv")
    if file_exists(read_stats_file):
        os.remove(read_stats_file)

    bam_file, seq_obj = _clean_alignment(args)

    logger.info("Parsing matrix file")
    seqL, y, l = parse_ma_file(seq_obj, args.ffile)
    # y, l = _total_counts(seqL.keys(), seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'aligned'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    if len(seqL.keys()) < 10:
        logger.error("It seems you have low coverage. Please check your fastq files have enough sequences.")
        raise ValueError("So few sequences.")

    logger.info("Cleaning bam file")
    y, l = _total_counts(seqL.keys(), seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'cleaned'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    clusL = _create_clusters(seqL, bam_file, args)
    y, l = _total_counts(clusL.seq.keys(), clusL.seq, aligned=True)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'clusters'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')

    logger.info("Solving multi-mapping events in the network of clusters")
    clusLred = _cleaning(clusL, args.dir_out)
    y, l = _total_counts(clusLred.clus, seqL)
    logger.info("counts after: %s" % sum(y.values()))
    logger.info("# sequences after: %s" % l)
    dt = pd.DataFrame({'sample': y.keys(), 'counts': y.values()})
    dt['step'] = 'meta-cluster'
    dt.to_csv(read_stats_file, sep="\t", index=False, header=False, mode='a')
    logger.info("Clusters up to %s" % (len(clusLred.clus.keys())))

    if args.show:
        logger.info("Creating sequences alignment to precursor")
        clusLred = show_seq(clusLred, args.index)

    clusLred = peak_calling(clusLred)

    clusLred = _annotate(args, clusLred)

    logger.info("Creating json and count matrix")
    json_file = _create_json(clusLred, args)

    logger.info("Output file in: %s" % args.dir_out)

    if args.db:
        name = args.db + ".db"
        logger.info("Create database: database/" + name)
        data = load_data(json_file)
        out_dir = op.join(args.dir_out, "database")
        make_database(data, name, out_dir)
    logger.info("Finished")
914010
def get_interested_subsequences(subsequences):#, subseqdim):
    
    keep_indices = []
    # if subseqdim == 'ALLD':
    #     filter_subseq_suffixes = filter_subseq_suffixes
    # elif subseqdim == '3D':
    #     filter_subseq_suffixes = filter_subseq_3D_suffixes
    # elif subseqdim == '2D':
    #     filter_subseq_suffixes = filter_subseq_2D_suffixes
    # else:
    #     raise ValueError('ALLD, 3D, or 2D only!')

    for subseq in subsequences:
        if subseq == 'all':
            keep_indices.extend([x for x in filter_subseq_suffixes[subseq]])
        else:
            keep_indices.extend([subseq + '_' + x for x in filter_subseq_suffixes[subseq]])
    return keep_indices
348044
def has_rabf_motif(self):
        

        if self.rabf_motifs:
            for gdomain in self.gdomain_regions:
                beg, end = map(int, gdomain.split('-'))
                motifs = [x for x in self.rabf_motifs if x[1] >= beg and x[2] <= end]
                if motifs:
                    matches = int(pairwise2.align.globalxx('12345', ''.join(str(x[0]) for x in motifs))[0][2])
                    if matches >= self.motif_number:
                        return True
        return False
63233
def highest(self):
        
        scores = self.scores()
        if not scores:
            return None
        maxscore = max(map(score, scores))
        return filter(lambda x: score(x) == maxscore, scores)
1173088
def _scan_smaller(self, seq, threshold=''):
        
        ll = self.ll #Shortcut for Log-likelihood matrix
        matches   = []
        endpoints = []
        scores    = []
        w         = self.width
        oseq      = seq
        seq       = seq.upper()
        for offset in range(self.width-len(seq)+1):    #Check if +/-1 needed
            maximum = 0
            for i in range(len(seq)):
                maximum = maximum + max(ll[i+offset].values())
            if not threshold: threshold = 0.8 * maximum
            total_f = 0
            total_r = 0
            for i in range(len(seq)):
                total_f = total_f + ll[i+offset      ][        seq[i] ]
                total_r = total_r + ll[w-(i+offset)-1][revcomp[seq[i]]]
            if 0:
                print "\t\t%s vs %s: F=%6.2f R=%6.2f %6.2f %4.2f"%(oseq, self.oneletter[offset:offset+len(seq)],
                                                                   total_f, total_r,  maximum,
                                                                   max([total_f,total_r])/self.maxscore)
            if total_f > threshold and total_f > total_r:
                endpoints.append( (offset,offset+self.width-1) )
                scores.append(total_f)
                matches.append(oseq[offset:offset+self.width])
            elif total_r > threshold:
                endpoints.append( (offset,offset+self.width-1) )
                scores.append(total_r)
                matches.append(oseq[offset:offset+self.width])
        return(matches,endpoints,scores)
318212
def findBest(self, pattern):
        
        findFailedRetry = True
        while findFailedRetry:
            best_match = None
            all_matches = self.findAll(pattern)
            for match in all_matches:
                if best_match is None or best_match.getScore() < match.getScore():
                    best_match = match
            self._lastMatch = best_match
            if best_match is not None:
                break
            path = pattern.path if isinstance(pattern, Pattern) else pattern
            findFailedRetry = self._raiseFindFailed("Could not find pattern '{}'".format(path))
            if findFailedRetry:
                time.sleep(self._repeatWaitTime)
        return best_match
616870
private int selectRightToLeft( int col , int[] scores ) {
		// see how far it can search
		int localMax = Math.min(imageWidth-regionWidth,col+maxDisparity)-col-minDisparity;

		int indexBest = 0;
		int indexScore = col;
		int scoreBest = scores[col];
		indexScore += imageWidth+1;

		for( int i = 1; i < localMax; i++ ,indexScore += imageWidth+1) {
			int s = scores[indexScore];

			if( s < scoreBest ) {
				scoreBest = s;
				indexBest = i;
			}
		}

		return indexBest;
	}
1335117
def hasScoreBetterThan(self, score):
        
        # Note: Do not assume that HSPs in an alignment are sorted in
        # decreasing order (as they are in BLAST output). If we could
        # assume that, we could just check the first HSP in each alignment.
        for hsp in self.hsps():
            if hsp.betterThan(score):
                return True
        return False
403260
def _kmp_search_all(self, pInput_sequence, pPattern):
                """"""
                r = []
                input_sequence, pattern = [len(bin(e)) for e in pInput_sequence], [len(bin(e)) for e in pPattern]
                n, m = len(input_sequence), len(pattern)
                d = p = 0
                next = self._kmp_construct_next(pattern)
                while d < n:
                        p = next[input_sequence[d]][p]
                        d += 1
                        if p == m:
                                r.append(d - m)
                                p = 0
                return r
664202
def makeSequenceAbsolute(relVSequence, minV, maxV):
    '''
    
    '''

    return [(value * (maxV - minV)) + minV for value in relVSequence]
648755
def top_matches(self, top):
        '''
        
        '''
        matches = DefaultOrderedDict(OrderedDict)
        # pylint: disable=cell-var-from-loop
        for saltenv, body in six.iteritems(top):
            if self.opts['saltenv']:
                if saltenv != self.opts['saltenv']:
                    continue
            for match, data in six.iteritems(body):
                def _filter_matches(_match, _data, _opts):
                    if isinstance(_data, six.string_types):
                        _data = [_data]
                    if self.matchers['confirm_top.confirm_top'](
                            _match,
                            _data,
                            _opts
                            ):
                        if saltenv not in matches:
                            matches[saltenv] = []
                        for item in _data:
                            if 'subfilter' in item:
                                _tmpdata = item.pop('subfilter')
                                for match, data in six.iteritems(_tmpdata):
                                    _filter_matches(match, data, _opts)
                            if isinstance(item, six.string_types):
                                matches[saltenv].append(item)
                            elif isinstance(item, dict):
                                env_key, inc_sls = item.popitem()
                                if env_key not in self.avail:
                                    continue
                                if env_key not in matches:
                                    matches[env_key] = []
                                matches[env_key].append(inc_sls)
                _filter_matches(match, data, self.opts['nodegroups'])
        ext_matches = self._master_tops()
        for saltenv in ext_matches:
            top_file_matches = matches.get(saltenv, [])
            if self.opts.get('master_tops_first'):
                first = ext_matches[saltenv]
                second = top_file_matches
            else:
                first = top_file_matches
                second = ext_matches[saltenv]
            matches[saltenv] = first + [x for x in second if x not in first]

        # pylint: enable=cell-var-from-loop
        return matches
635749
def fqscreen_simple_plot(self):
        

        # First, sum the different types of alignment counts
        data = OrderedDict()
        cats = OrderedDict()
        for s_name in sorted(self.fq_screen_data):
            data[s_name] = OrderedDict()
            sum_alignments = 0
            for org in self.fq_screen_data[s_name]:
                if org == 'total_reads':
                    continue
                try:
                    data[s_name][org] = self.fq_screen_data[s_name][org]['counts']['one_hit_one_library']
                except KeyError:
                    log.error("No counts found for '{}' ('{}'). Could be malformed or very old FastQ Screen results.".format(org, s_name))
                    continue
                try:
                    data[s_name][org] += self.fq_screen_data[s_name][org]['counts']['multiple_hits_one_library']
                except KeyError:
                    pass
                sum_alignments += data[s_name][org]
                if org not in cats and org != 'No hits':
                    cats[org] = { 'name': org }

            # Calculate hits in multiple genomes
            if 'total_reads' in self.fq_screen_data[s_name]:
                data[s_name]['Multiple Genomes'] = self.fq_screen_data[s_name]['total_reads'] - sum_alignments

        # Strip empty dicts
        for s_name in list(data.keys()):
            if len(data[s_name]) == 0:
                del data[s_name]

        pconfig = {
            'id': 'fastq_screen',
            'title': 'FastQ Screen',
            'cpswitch_c_active': False,
            'hide_zero_cats': False
        }
        cats['Multiple Genomes'] = { 'name': 'Multiple Genomes', 'color': '#820000' }
        cats['No hits'] = { 'name': 'No hits', 'color': '#cccccc' }

        return bargraph.plot(data, cats, pconfig)
774982
def make_simple_col_consensus(bg_freqs):
    
    # Hack: use default kwargs to persist across iterations
    def col_consensus(col, prev_col=[], prev_char=[]):
        # Count the amino acid types in this column
        aa_counts = sequtils.aa_frequencies(col)
        assert aa_counts, "Column is all gaps! That's not allowed."
        # Take the most common residue(s)
        best_char, best_score = max(aa_counts.iteritems(),
                                    key=lambda kv: kv[1])
        # Resolve ties
        ties = [aa for aa in aa_counts if aa_counts[aa] == best_score]
        if len(ties) > 1:
            # Breaker #1: most common after the prev. consensus char
            # Resolve a tied col by restricting to rows where the preceding
            # char is the consensus type for that (preceding) col
            if prev_char and prev_col:
                mc_next = Counter(
                        [b for a, b in zip(prev_col, col)
                            if a == prev_char[0] and b in ties]
                        ).most_common()
                ties_next = [x[0] for x in mc_next
                        if x[1] == mc_next[0][1]]
                if ties_next:
                    ties = ties_next
            if len(ties) > 1:
                # Breaker #2: lowest overall residue frequency
                ties.sort(key=lambda aa: bg_freqs[aa])
            best_char = ties[0]
        else:
            assert best_char == ties[0], \
                    'WTF %s != %s[0]' % (best_char, ties)
        # Save values for tie-breaker #1
        prev_col[:] = col
        prev_char[:] = best_char
        return best_char
    return col_consensus
1212960
def match_msequence(self, tokens, item):
        """"""
        series_type, head_matches, middle, _, last_matches = tokens
        self.add_check("_coconut.isinstance(" + item + ", _coconut.abc.Sequence)")
        self.add_check("_coconut.len(" + item + ") >= " + str(len(head_matches) + len(last_matches)))
        if middle != wildcard:
            if len(head_matches) and len(last_matches):
                splice = "[" + str(len(head_matches)) + ":" + str(-len(last_matches)) + "]"
            elif len(head_matches):
                splice = "[" + str(len(head_matches)) + ":]"
            elif len(last_matches):
                splice = "[:" + str(-len(last_matches)) + "]"
            else:
                splice = ""
            self.assign_to_series(middle, series_type, item + splice)
        self.match_all_in(head_matches, item)
        for i, match in enumerate(last_matches):
            self.match(match, item + "[" + str(i - len(last_matches)) + "]")
1346834
def best_match_from_list(item,options,fuzzy=90,fname_match=True,fuzzy_fragment=None,guess=False):
    ''''''
    matches = matches_from_list(item,options,fuzzy,fname_match,fuzzy_fragment,guess)
    if len(matches)>0:
        return matches[0]
    return None
1219989
def rank_motifs(stats, metrics=("roc_auc", "recall_at_fdr")):
    """"""
    rank = {}
    combined_metrics = []
    motif_ids = stats.keys()
    background = list(stats.values())[0].keys()
    for metric in metrics:
        mean_metric_stats = [np.mean(
            [stats[m][bg][metric] for bg in background]) for m in motif_ids]
        ranked_metric_stats = rankdata(mean_metric_stats)
        combined_metrics.append(ranked_metric_stats)
    
    for motif, val in zip(motif_ids, np.mean(combined_metrics, 0)):
        rank[motif] = val

    return rank