def best_kmers(dt, response, sequence, k=6, consider_shift=True, n_cores=1,
               seq_align="start", trim_seq_len=None):
    y = dt[response]
    seq = dt[sequence]
    if trim_seq_len is not None:
        seq = pad_sequences(seq, align=seq_align, maxlen=trim_seq_len)
        seq = [s.replace("N", "") for s in seq]
    dt_kmer = kmer_count(seq, k)
    Xsp = csc_matrix(dt_kmer)
    en = ElasticNet(alpha=1, standardize=False, n_splits=3)
    en.fit(Xsp, y)
    nonzero_kmers = dt_kmer.columns.values[en.coef_ != 0].tolist()
    def find_next_best(dt_kmer, y, selected_kmers, to_be_selected_kmers, consider_shift=True):
        F, pval = f_regression(dt_kmer[to_be_selected_kmers], y)
        kmer = to_be_selected_kmers.pop(pval.argmin())
        selected_kmers.append(kmer)
        def select_criterion(s1, s2, consider_shift=True):
            if hamming_distance(s1, s2) <= 1:
                return False
            if consider_shift and hamming_distance(s1[1:], s2[:-1]) == 0:
                return False
            if consider_shift and hamming_distance(s1[:-1], s2[1:]) == 0:
                return False
            return True
        to_be_selected_kmers = [ckmer for ckmer in to_be_selected_kmers
                                if select_criterion(ckmer, kmer, consider_shift)]
        if len(to_be_selected_kmers) == 0:
            return selected_kmers
        else:
            lm = LinearRegression()
            lm.fit(dt_kmer[selected_kmers], y)
            y_new = y - lm.predict(dt_kmer[selected_kmers])
            return find_next_best(dt_kmer, y_new, selected_kmers, to_be_selected_kmers, consider_shift)
    selected_kmers = find_next_best(dt_kmer, y, [], nonzero_kmers, consider_shift)
    return selected_kmers