La tarea es contar el número de subcadenas distintas de longitud k, para k = 1,2,3,4, .....

Salida

Debe generar una línea por k logro que complete con un número por línea de salida. Su producción debe estar en orden de aumento khasta que se acabe el tiempo.

Puntuación

Su puntaje es el k más alto que puede alcanzar en mi computadora en menos de 1 minuto.

Debe usar http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/chr2.fa.gz como su entrada e ignorar las nuevas líneas. Sin embargo, su código debe ser sensible a mayúsculas y minúsculas.

Puede descomprimir la entrada antes de comenzar la sincronización.

El siguiente código (ineficiente) cuenta el número de 4 mers distintos.

awk -v substr_length=4 '(len=length($0))>=substr_length{for (i=1; (i-substr_length)<len; i++) substrs[substr($0,i,substr_length)]++}; END{for (i in substrs) print substrs[i], i}' file.txt|wc

Límites de memoria

Para hacer que su código se comporte bien en mi computadora y para que la tarea sea más desafiante, limitaré la RAM que puede usar a 2GB usando

ulimit -v 2000000

antes de ejecutar su código. Soy consciente de que esta no es una forma sólida de limitar el uso de RAM, así que no utilice formas imaginativas para superar este límite, por ejemplo, generando nuevos procesos. Por supuesto, puede escribir código multiproceso, pero si alguien lo hace, tendré que aprender a limitar la RAM total utilizada para eso.

Desempate

En el caso de un empate por un máximo k, calcularé cuánto tiempo se tarda en dar los resultados k+1y gana el más rápido. En el caso de que se ejecuten al mismo tiempo hasta dentro de un segundo k+1, la primera presentación gana.

Idiomas y bibliotecas

Puede usar cualquier idioma que tenga un compilador / intérprete / etc. para Linux y cualquier biblioteca que también esté disponible gratuitamente para Linux.

Mi máquina Los tiempos se ejecutarán en mi máquina. Esta es una instalación estándar de ubuntu en un procesador AMD FX-8350 de ocho núcleos en una placa base Asus M5A78L-M / USB3 (Socket AM3 +, 8GB DDR3). Esto también significa que necesito poder ejecutar su código. Como consecuencia, solo use software gratuito fácilmente disponible e incluya instrucciones completas sobre cómo compilar y ejecutar su código.

Prueba de salida

El código de FUZxxl genera lo siguiente (pero no todo en 1 minuto) que creo que es correcto.

14
92
520
2923
15714
71330
265861
890895
2482912
5509765
12324706
29759234

tabla de la liga

• k> = 4000 FUZxxl (C)
• k = 16 por Keith Randall (C ++)
• k = 10 por FUZxxl (C)

¿Cuánto puede especializar su código a la entrada?

• Claramente, arruinaría la competencia si solo calculara previamente las respuestas y su código las enviara. No hagas eso.
• Idealmente, cualquier cosa que su código necesite para aprender sobre los datos que usará para ejecutarse más rápidamente, puede aprender en tiempo de ejecución.
• Sin embargo, puede suponer que la entrada se verá como los datos en los archivos * .fa en http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/ .

C (≥ 4000)

Este código no usa menos de 2 GB de RAM (usa un poco más) ni produce ningún resultado en el primer minuto. Pero si espera unos seis minutos, imprime todos los recuentos de k -mer a la vez.

Se incluye una opción para limitar la k más alta para la cual contamos los k -mers. Cuando k se limita a un rango razonable, el código termina en menos de un minuto y usa menos de 2 GiB de RAM. OP ha calificado esta solución y finaliza en menos de un minuto para un límite no significativamente superior a 4000.

¿Como funciona?

El algoritmo tiene cuatro pasos.

1. Lea el archivo de entrada en un búfer y elimine las nuevas líneas.

2. Sufijo-ordena una matriz de índices en el búfer de entrada. Por ejemplo, los sufijos de la cadena mississippison:

   mississippi
   ississippi
   ssissippi
   sissippi
   issippi
   ssippi
   sippi
   ippi
   ppi
   pi
   i
   

   Estas cadenas ordenadas en orden lexicográfico son:

   i
   ippi
   issippi
   ississippi
   mississippi
   pi
   ppi
   sippi
   sissippi
   ssippi
   ssissippi
   

   Es fácil ver que todas las subcadenas iguales de longitud k para todas las k se encuentran en entradas adyacentes de la matriz ordenada por sufijo.

3. Se rellena una matriz entera en la que almacenamos en cada índice k el número de k -mers distintos . Esto se hace de una manera ligeramente complicada para acelerar el proceso. Considere dos entradas adyacentes como la matriz ordenada por sufijo.

      p     l
      v     v
   issippi
   ississippi
   

   p denota la longitud del prefijo común más largo de las dos entradas, l denota la longitud de la segunda entrada. Para tal par, encontramos una nueva subcadena distinta de longitud k para p < k ≤ l . Como p ≪ l a menudo se mantiene, no es práctico incrementar una gran cantidad de entradas de matriz para cada par. En su lugar, almacenamos la matriz como una matriz de diferencia, donde cada entrada k denota la diferencia al número de k -mers y el número de ( k  - 1) -mers. Esto convierte una actualización del formulario

   0  0  0  0 +1 +1 +1 +1 +1 +1  0  0  0
   

   en una actualización mucho más rápida del formulario

   0  0  0  0 +1  0  0  0  0  0 -1  0  0
   

   Al observar cuidadosamente que l siempre es diferente y, de hecho, cada 0 < l < n aparecerá exactamente una vez, podemos omitir las restas y restar 1 de cada diferencia al convertir las diferencias en cantidades.

4. Las diferencias se convierten en importes y se imprimen.

El código fuente

Esta fuente utiliza libdivsufsort para ordenar las matrices de sufijos. El código genera resultados de acuerdo con la especificación cuando se compila con esta invocación.

cc -O -o dsskmer dsskmer.c -ldivsufsort

alternativamente, el código puede generar una salida binaria cuando se compila con la siguiente invocación.

cc -O -o dsskmer -DBINOUTPUT dsskmer.c -ldivsufsort

Para limitar la k más alta para la cual se deben contar k -mers, suministre -DICAP = k donde k es el límite:

cc -O -o dsskmer -DICAP=64 dsskmer.c -ldivsufsort

Compile con -O3si su compilador proporciona esta opción.

#include <stdlib.h>
#include <stdio.h>

#include <divsufsort.h>

#ifndef BINOUTPUT
static char stdoutbuf[1024*1024];
#endif

/*
 * load input from stdin and remove newlines. Save length in flen.
 */
static unsigned char*
input(flen)
    int *flen;
{
    size_t n, len, pos, i, j;
    off_t slen;
    unsigned char *buf, *sbuf;

    if (fseek(stdin, 0L, SEEK_END) != 0) {
        perror("Cannot seek stdin");
        abort();
    }

    slen = ftello(stdin);
    if (slen == -1) {
        perror("Cannot tell stdin");
        abort();
    }

    len = (size_t)slen;

    rewind(stdin);

    /* Prepare for one extra trailing \0 byte */
    buf = malloc(len + 1);
    if (buf == NULL) {
        perror("Cannot malloc");
        abort();
    }

    pos = 0;

    while ((n = fread(buf + pos, 1, len - pos, stdin)) != 0)
        pos += n;

    if (ferror(stdin)) {
        perror("Cannot read from stdin");
        abort();
    }

    /* remove newlines */
    for (i = j = 0; i < len; i++)
        if (buf[i] != '\n')
            buf[j++] = buf[i];

    /* try to reclaim some memory */
    sbuf = realloc(buf, j);
    if (sbuf == NULL)
        sbuf = buf;

    *flen = (int)j;
    return sbuf;
}

/*
 * Compute for all k the number of k-mers. kmers will contain at index i the
 * number of (i + 1) mers. The count is computed as an array of differences,
 * where kmers[i] == kmersum[i] - kmersum[i-1] + 1 and then summed up by the
 * caller. This algorithm is a little bit unclear, but when you write subsequent
 * suffixes of an array on a piece of paper, it's easy to see how and why it
 * works.
 */
static void
count(buf, sa, kmers, n)
    const unsigned char *buf;
    const int *sa;
    int *kmers;
{
    int i, cl, cp;

    /* the first item needs special treatment */
    kmers[0]++;

    for (i = 1; i < n; i++) {
        /* The longest common prefix of the two suffixes */
        cl = n - (sa[i-1] > sa[i] ? sa[i-1] : sa[i]);
#ifdef ICAP
        cl = (cl > ICAP ? ICAP : cl);
#endif
        for (cp = 0; cp < cl; cp++)
            if (buf[sa[i-1] + cp] != buf[sa[i] + cp])
                break;

        /* add new prefixes to the table */
        kmers[cp]++;
    }
}

extern int
main()
{
    unsigned char *buf;
    int blen, ilen, *sa, *kmers, i;

    buf = input(&blen);

    sa = malloc(blen * sizeof *sa);

    if (divsufsort(buf, sa, blen) != 0) {
        puts("Cannot divsufsort");
        abort();
    }

#ifdef ICAP
    ilen = ICAP;
    kmers = calloc(ilen + 1, sizeof *kmers);
#else
    ilen = blen;
    kmers = calloc(ilen, sizeof *kmers);
#endif

    if (kmers == NULL) {
        perror("Cannot malloc");
        abort();
    }

    count(buf, sa, kmers, blen);

#ifndef BINOUTPUT
    /* sum up kmers differences */
    for (i = 1; i < ilen; i++)
        kmers[i] += kmers[i-1] - 1;


    /* enlarge buffer of stdout for better IO performance */
    setvbuf(stdout, stdoutbuf, _IOFBF, sizeof stdoutbuf);

    /* human output */
    for (i = 0; i < ilen; i++)
        printf("%d\n", kmers[i]);
#else
    /* binary output in host endianess */
    fprintf(stderr, "writing out result...\n");
    fwrite(kmers, sizeof *kmers, ilen, stdout);
#endif

    return 0;
}

El formato de archivo binario se puede convertir al formato de salida legible por humanos con el siguiente programa:

#include <stdlib.h>
#include <stdio.h>

static int inbuf[BUFSIZ];
static char outbuf[BUFSIZ];

extern int main()
{
    int i, n, sum = 1;

    setbuf(stdout, outbuf);

    while ((n = fread(inbuf, sizeof *inbuf, BUFSIZ, stdin)) > 0)
        for (i = 0; i < n; i++)
            printf("%d\n", sum += inbuf[i] - 1);

    if (ferror(stdin)) {
        perror("Error reading input");
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}

salida de muestra

Aquíchr22.fa puede encontrar una salida de muestra en formato binario para el archivo . Descomprima con primero. La salida se proporciona en formato binario solo porque se comprime mucho mejor (3,5 kB frente a 260 MB) que la salida en formato legible para humanos. Sin embargo, tenga en cuenta que la salida de referencia tiene un tamaño de 924 MB cuando se descomprime. Es posible que desee utilizar una tubería como esta:bzip2 -d

bzip2 -dc chr2.kmerdiffdat.bz2 | ./unbin | less

C ++, k = 16, 37 segundos

Calcula todos los 16 mers en la entrada. Cada 16 meros se empaqueta 4 bits a un símbolo en una palabra de 64 bits (con un patrón de un bit reservado para EOF). Las palabras de 64 bits se ordenan. La respuesta para cada k puede leerse observando con qué frecuencia cambian los 4 * k bits superiores de las palabras ordenadas.

Para la entrada de prueba, uso aproximadamente 1.8 GB, justo debajo del cable.

Estoy seguro de que la velocidad de lectura podría mejorarse.

Salida:

14
92
520
2923
15714
71330
265861
890895
2482912
5509765
12324706
29759234
69001539
123930801
166196504
188354964

Programa:

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>

int main(int argc, char *argv[]) {
    // read file
    printf("reading\n");
    FILE *f = fopen(argv[1], "rb");
    fseek(f, 0, SEEK_END);
    int32_t size = ftell(f);
    printf("size: %d\n", size);
    fseek(f, 0, SEEK_SET);

    // table to convert 8-bit input into 4-bit tokens.  Reserve 15 for EOF.
    int ntokens = 0;
    int8_t squash[256];
    for(int i = 0; i < 256; i++) squash[i] = -1;

    uint64_t *buf = (uint64_t*)malloc(8*size);

    int32_t n = 0;
    uint64_t z = 0;
    for(int32_t i = 0; i < size; i++) {
        char c = fgetc(f);
        if(c == '\n') continue;
        int8_t s = squash[c];
        if(s == -1) {
            if(ntokens == 15) {
                printf("too many tokens\n");
                exit(1);
            }
            squash[c] = s = ntokens++;
        }
        z <<= 4;
        z += s;
        n++;

        if(n >= 16) buf[n-16] = z;
    }
    for(int32_t i = 1; i < 16; i++) {
        z <<= 4;
        z += 15;
        buf[n-16+i] = z;
    }
    printf("   n: %d\n", n);

    // sort these uint64_t's
    printf("sorting\n");
    std::sort(buf, buf+n);

    for(int32_t k = 1; k <= 16; k++) {
        // count unique entries
        int32_t shift = 64-4*k;
        int32_t cnt = 1;
        int64_t e = buf[0] >> shift;
        for(int32_t i = 1; i < n; i++) {
            int64_t v = buf[i] >> shift;
            if((v & 15) == 15) continue; // ignore EOF entries
            if(v != e) {
                cnt++;
                e = v;
            }
        }

        printf("%d\n", cnt);
    }
}

Compilar g++ -O3 kmer.cc -o kmery ejecutar con ./kmer chr2.fa.

C ++: una mejora con respecto a la solución FUZxxl

No merezco absolutamente ningún crédito por el método de cálculo en sí, y si no se presenta un mejor enfoque mientras tanto, la recompensa debería ir a FUZxxl por derecho.

#define _CRT_SECURE_NO_WARNINGS // a Microsoft thing about strcpy security issues
#include <vector>

#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cstring>
using namespace std;

#include "divsufsort.h"

// graceful exit of sorts
void panic(const char * msg)
{
    cerr << msg;
    exit(0);
}

// approximative timing of various steps
struct tTimer {
    time_t begin;
    tTimer() { begin = time(NULL); }
    void print(const char * msg)
    {
        time_t now = time(NULL);
        cerr << msg << " in " << now - begin << "s\n";
        begin = now;
    }
};

// load input pattern
unsigned char * read_sequence (const char * filename, int& len)
{
    ifstream file(filename);
    if (!file) panic("could not open file");

    string str;
    std::string line;
    while (getline(file, line)) str += line;
    unsigned char * res = new unsigned char[str.length() + 1];
    len = str.length()+1;
    strcpy((char *)res, str.c_str());
    return res;
}

#ifdef FUZXXL_METHOD
/*
* Compute for all k the number of k-mers. kmers will contain at index i the
* number of (i + 1) mers. The count is computed as an array of differences,
* where kmers[i] == kmersum[i] - kmersum[i-1] + 1 and then summed up by the
* caller. This algorithm is a little bit unclear, but when you write subsequent
* suffixes of an array on a piece of paper, it's easy to see how and why it
* works.
*/
static void count(const unsigned char *buf, const int *sa, int *kmers, int n)
{
    int i, cl, cp;

    /* the first item needs special treatment */
    /*
        kuroi neko: since SA now includes the null string, kmers[0] is indeed 0 instead of 1
    */
//  kmers[0]++;

    for (i = 1; i < n; i++) {
        /* The longest common prefix of the two suffixes */
        cl = n - (sa[i - 1] > sa[i] ? sa[i - 1] : sa[i]);
#ifdef ICAP
        cl = (cl > ICAP ? ICAP : cl);
#endif
        for (cp = 0; cp < cl; cp++)
        if (buf[sa[i - 1] + cp] != buf[sa[i] + cp])
            break;

        /* add new prefixes to the table */
        kmers[cp]++;
    }
}

#else // Kasai et al. method

// compute kmer cumulative count using Kasai et al. LCP construction algorithm
void compute_kmer_cumulative_sums(const unsigned char * t, const int * sa, int * kmer, int len)
{
    // build inverse suffix array
    int * isa = new int[len];
    for (int i = 0; i != len; i++) isa[sa[i]] = i;

    // enumerate common prefix lengths
    memset(kmer, 0, len*sizeof(*kmer));
    int lcp = 0;
    int limit = len - 1;
    for (int i = 0; i != limit; i++)
    {
        int k = isa[i];
        int j = sa[k - 1];
        while (t[i + lcp] == t[j + lcp]) lcp++;

        // lcp now holds the kth longest commpn prefix length, which is just what we need to compute our kmer counts
        kmer[lcp]++;
        if (lcp > 0) lcp--;
    }
    delete[] isa;
}

#endif // FUZXXL_METHOD

int main (int argc, char * argv[])
{
    if (argc != 2) panic ("missing data file name");

    tTimer timer;
    int blen;
    unsigned char * sequence;
    sequence = read_sequence(argv[1], blen);
    timer.print("input read");

    vector<int>sa;
    sa.assign(blen, 0);
    if (divsufsort(sequence, &sa[0], blen) != 0) panic("divsufsort failed");
    timer.print("suffix table constructed");

    vector<int>kmers;
    kmers.assign(blen,0);

#ifdef FUZXXL_METHOD
    count(sequence, &sa[0], &kmers[0], blen);
    timer.print("FUZxxl count done");
#else
    compute_kmer_cumulative_sums(sequence, &sa[0], &kmers[0], blen);
    timer.print("Kasai  count done");
#endif

    /* sum up kmers differences */
    for (int i = 1; i < blen; i++) kmers[i] += kmers[i - 1] - 1;
    timer.print("sum done");

    /* human output */

    if (blen>10) blen = 10; // output limited to the first few values to avoid cluttering display or saturating disks

    for (int i = 0; i != blen; i++) printf("%d ", kmers[i]);
    return 0;
}

Simplemente usé Kasai et al. algoritmo para calcular LCP en O (n).
El resto es una mera adaptación del código FUZxxl, utilizando características C ++ más concisas aquí y allá.

Dejé el código de cálculo original para permitir comparaciones.

Como los procesos más lentos son la construcción de SA y el recuento de LCP, eliminé la mayoría de las otras optimizaciones para evitar saturar el código para obtener ganancias despreciables.

Extendí la tabla SA para incluir el prefijo de longitud cero. Eso facilita el cálculo de LCP.

No proporcioné una opción de limitación de longitud, el proceso más lento ahora es el cálculo de SA que no puede reducirse (o al menos no veo cómo podría ser).

También eliminé la opción de salida binaria y limité la visualización a los primeros 10 valores.
Supongo que este código es solo una prueba de concepto, por lo que no es necesario saturar las pantallas o saturar los discos.

Construyendo el ejecutable

Tuve que compilar todo el proyecto (incluida la versión lite dedivsufsort ) para x64 para superar el límite de asignación de Win32 de 2 Gb.

divsufsortel código arroja un montón de advertencias debido al uso intensivo de ints en lugar de size_ts, pero eso no será un problema para las entradas de menos de 2 Gb (que de todos modos requerirían 26 Gb de RAM: D).

Construcción de Linux

compilar main.cppy divsufsort.cusar los comandos:

g++ -c -O3 -fomit-frame-pointer divsufsort.c 
g++ -O3 main.cpp -o kmers divsufsort.o

Supongo que la divsufsortbiblioteca normal debería funcionar bien en Linux nativo, siempre y cuando pueda asignar un poco más de 3Gb.

Actuaciones

El algoritmo Kasai requiere la tabla SA inversa, que consume 4 bytes más por carácter para un total de 13 (en lugar de 9 con el método FUZxxl).

El consumo de memoria para la entrada de referencia es, por lo tanto, superior a 3Gb.

Por otro lado, el tiempo de cálculo se mejora dramáticamente, y todo el algoritmo ahora está en O (n):

input read in 5s
suffix table constructed in 37s
FUZxxl count done in 389s
Kasai  count done in 27s
14 92 520 2923 15714 71330 265861 890895 2482912 5509765 (etc.)

Futuras mejoras

La construcción de SA es ahora el proceso más lento.
Algunos bits del divsufsortalgoritmo están destinados a ser paralelos con cualquier característica incorporada de un compilador desconocido para mí, pero si es necesario, el código debería ser fácil de adaptar a subprocesos múltiples más clásicos ( por ejemplo, a la C ++ 11).
La biblioteca también tiene una gran cantidad de parámetros, incluidos varios tamaños de cubetas y la cantidad de símbolos distintos en la cadena de entrada. Solo les eché un vistazo superficial, pero sospecho que valdría la pena intentar comprimir el alfabeto si tus cadenas son interminables permutaciones de ACTG ( y estás desesperado por las actuaciones).

Existen algunos métodos paralelizables para calcular LCP desde SA también, pero dado que el código debería ejecutarse en menos de un minuto en un procesador un poco más potente que mi pequeño [email protected] y todo el algoritmo está en O (n), dudo que esto Valdría la pena el esfuerzo.

C (puede resolver hasta 10 en un minuto en mi máquina)

Esta es una solución muy simple. Construye un árbol de los k -mers encontrados y los cuenta. Para conservar la memoria, los caracteres se convierten primero en enteros de 0 a n - 1 donde n es el número de caracteres diferentes en la entrada, por lo que no necesitamos proporcionar espacio para los caracteres que nunca aparecen. Además, se asigna menos memoria para las hojas que para otros nodos para conservar más memoria. Esta solución utiliza alrededor de 200 MiB de RAM durante su tiempo de ejecución en mi máquina. Todavía lo estoy mejorando, ¡así que quizás en la iteración podría ser aún más rápido!

Para compilar, guarde el siguiente código en un archivo llamado kmers.cy luego ejecútelo en un sistema operativo similar a POSIX:

cc -O -o kmers kmers.c

Es posible que desee sustituir -O3por -Osi sus soportes compilador eso. Para ejecutar, desembalar por primera vez chr2.fa.gzen chr2.fay vuelva a ejecutar:

./kmers <chr2.fa

Esto produce la salida paso a paso. Es posible que desee limitar tanto el tiempo como el espacio. Usa algo como

( ulimit -t 60 -v 2000000 ; ./kmers <chrs.fa )

para reducir recursos según sea necesario.

#include <limits.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

/*
 * A tree for the k-mers. It contains k layers where each layer is an
 * array of struct tr. The last layer contains a 1 for every kmer found,
 * the others pointers to the next layer or NULL.
 */
struct tr {
    struct tr *n;
};

/* a struct tr to point to */
static struct tr token;

static void     *mem(size_t s);
static int       add(struct tr*, unsigned char*, int, size_t);
static unsigned char    *input(size_t*);
extern int       main(void);

/*
 * Allocate memory, fail if we can't.
 */
static void*
mem(s)
    size_t s;
{
    void *p;

    p = calloc(s, 1);
    if (p != NULL)
        return p;

    perror("Cannot malloc");
    abort();
}

/*
 * add s to b, return 1 if added, 0 if present before. Assume that n - 1 layers
 * of struct tr already exist and only add an n-th layer if needed. In the n-th
 * layer, a NULL pointer denotes non-existance, while a pointer to token denotes
 * existance.
 */
static int
add(b, s, n, is)
    struct tr *b;
    unsigned char *s;
    size_t is;
{
    struct tr **nb;
    int added;
    int i;

    for (i = 0; i < n - 2; i++) {
        b = b[s[i]].n;
    }

    nb = &b[s[n - 2]].n;

    if (*nb == NULL || *nb == &token)
        *nb = mem(is * sizeof *b);

    added = (*nb)[s[n - 1]].n == NULL;
    (*nb)[s[n - 1]].n = &token;

    return (added);
}

/*
 * load input from stdin and remove newlines. Save length in flen.
 */
static unsigned char*
input(flen)
    size_t *flen;
{
    size_t n, len, pos, i, j;
    unsigned char *buf;

    if (fseek(stdin, 0L, SEEK_END) != 0) {
        perror("Cannot seek stdin");
        abort();
    }

    len = ftello(stdin);
    if (len == -1) {
        perror("Cannot tell stdin");
        abort();
    }

    rewind(stdin);

    /* no need to zero out, so no mem() */
    buf = malloc(len);
    if (buf == NULL) {
        perror("Cannot malloc");
        abort();
    }

    pos = 0;

    while ((n = fread(buf + pos, 1, len - pos, stdin)) != 0)
        pos += n;

    if (ferror(stdin)) {
        perror("Cannot read from stdin");
        abort();
    }

    /* remove newlines */
    for (i = j = 0; i < len; i++)
        if (buf[i] != '\n')
            buf[j++] = buf[i];

    *flen = j;
    return buf;
}

extern int
main()
{
    struct tr *b;
    size_t flen, c, i, k, is;
    unsigned char *buf, itab[1 << CHAR_BIT];

    buf = input(&flen);

    memset(itab, 0, sizeof itab);

    /* process 1-mers */
    for (i = 0; i < flen; i++)
        itab[buf[i]] = 1;

    is = 0;
    for (i = 0; i < sizeof itab / sizeof *itab; i++)
        if (itab[i] != 0)
            itab[i] = is++;

    printf("%zd\n", is);

    /* translate characters into indices */
    for (i = 0; i < flen; i++)
        buf[i] = itab[buf[i]];

    b = mem(is * sizeof *b);

    /* process remaining k-mers */
    for (k = 2; k < flen; k++) {
        c = 0;

        for (i = 0; i < flen - k + 1; i++)
            c += add(b, buf + i, k, is);

        printf("%zd\n", c);
    }

    return 0;
}

Mejoras

1. 8 → 9: Lea todo el archivo al principio, preprocese una vez y manténgalo en el núcleo. Esto mejora enormemente el rendimiento.
2. Use menos memoria, escriba la salida correcta.
3. Repara el formato de salida nuevamente.
4. Arreglar off-by-one.
5. 9 → 10: No tires lo que ya hiciste.