Die Aufgabe besteht darin, die Anzahl der unterschiedlichen Teilzeichenfolgen der Länge k für k = 1,2,3,4, ..... zu zählen.

Ausgabe

Sie sollten eine Zeile pro Zeile ausgeben, die kSie mit einer Nummer pro Ausgabezeile vervollständigen können. Ihre Ausgabe sollte in der Reihenfolge steigen, kbis Ihnen die Zeit ausgeht.

Ergebnis

Ihre Punktzahl ist die höchste k, die Sie in weniger als 1 Minute auf meinem Computer erreichen können.

Sie sollten http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/chr2.fa.gz als Eingabe verwenden und Zeilenumbrüche ignorieren. Ihr Code sollte jedoch zwischen Groß- und Kleinschreibung unterscheiden.

Sie können den Eingang dekomprimieren, bevor Sie das Timing starten.

Der folgende (ineffiziente) Code zählt die Anzahl der verschiedenen 4-mers.

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

Speichergrenzen

Damit sich Ihr Code auf meinem Computer gut verhält und die Aufgabe schwieriger wird, beschränke ich den Arbeitsspeicher, den Sie verwenden können, auf 2 GB

ulimit -v 2000000

bevor Sie Ihren Code ausführen. Ich bin mir bewusst, dass dies keine solide Methode ist, um die RAM-Nutzung einzuschränken. Verwenden Sie daher bitte keine einfallsreichen Methoden, um diese Grenze zu umgehen, indem Sie beispielsweise neue Prozesse erzeugen. Natürlich können Sie Multithread-Code schreiben, aber wenn jemand dies tut, muss ich lernen, wie man den dafür verwendeten RAM begrenzt.

Tie Breaker

Im Falle eines Unentschieden für ein Maximum werde kich die Zeit bestimmen, wie lange es dauert, bis die Ergebnisse erreicht sind k+1und der schnellste gewinnt. Für den Fall, dass sie innerhalb einer Sekunde bis zur gleichen Zeit laufen k+1, gewinnt die erste Einreichung.

Sprachen und Bibliotheken

Sie können jede Sprache verwenden, die einen frei verfügbaren Compiler / Interpreter / etc. für Linux und alle Bibliotheken, die auch für Linux frei verfügbar sind.

Meine Maschine Die Timings werden auf meiner Maschine ausgeführt. Dies ist eine Standard-Ubuntu-Installation auf einem AMD FX-8350 Eight-Core-Prozessor auf einem Asus M5A78L-M / USB3-Motherboard (Sockel AM3 +, 8 GB DDR3). Dies bedeutet auch, dass ich Ihren Code ausführen kann. Verwenden Sie daher nur leicht verfügbare freie Software und fügen Sie bitte vollständige Anweisungen zum Kompilieren und Ausführen Ihres Codes bei.

Testausgabe

Der Code von FUZxxl gibt Folgendes aus (aber nicht alle innerhalb von 1 Minute), was ich für richtig halte.

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

Ligatabelle

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

Wie stark können Sie Ihren Code auf die Eingabe spezialisieren?

• Es wäre klar, dass dies die Konkurrenz ruinieren würde, wenn Sie nur die Antworten vorberechnen und von Ihrem Code ausgeben lassen würden. Tu das nicht.
• Im Idealfall muss alles, was Ihr Code benötigt, um mehr über die Daten zu erfahren, die zur schnelleren Ausführung verwendet werden, zur Laufzeit gelernt werden.
• Sie können jedoch davon ausgehen, dass die Eingabe wie die Daten in den * .fa-Dateien unter http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/ aussieht .

C (≥ 4000)

Dieser Code verwendet weder weniger als 2 GB RAM (er verbraucht etwas mehr) noch erzeugt er in der ersten Minute eine Ausgabe. Aber wenn Sie für etwa sechs Minuten warten, es druckt alle die k - mer zählt auf einmal.

Eine Option ist enthalten, um das höchste k zu begrenzen, für das wir die k- mere zählen. Wenn k auf einen vernünftigen Bereich begrenzt ist, endet der Code in weniger als einer Minute und verwendet weniger als 2 GiB RAM. OP hat diese Lösung bewertet und endet in weniger als einer Minute für einen Grenzwert, der nicht wesentlich höher als 4000 ist.

Wie funktioniert es?

Der Algorithmus besteht aus vier Schritten.

1. Lesen Sie die Eingabedatei in einen Puffer und entfernen Sie Zeilenumbrüche.

2. Sortieren Sie ein Array von Indizes nach dem Eingabepuffer. Die Suffixe der Zeichenfolge mississippilauten beispielsweise:

   mississippi
   ississippi
   ssissippi
   sissippi
   issippi
   ssippi
   sippi
   ippi
   ppi
   pi
   i
   

   Diese in lexikografischer Reihenfolge sortierten Zeichenfolgen sind:

   i
   ippi
   issippi
   ississippi
   mississippi
   pi
   ppi
   sippi
   sissippi
   ssippi
   ssissippi
   

   Es ist leicht zu erkennen, dass alle gleichen Teilzeichenfolgen der Länge k für alle k in benachbarten Einträgen des nach Suffix sortierten Arrays gefunden werden.

3. Es wird ein ganzzahliges Array ausgefüllt, in dem wir an jedem Index k die Anzahl der verschiedenen k- mere speichern. Dies geschieht auf leicht gewundene Weise, um den Prozess zu beschleunigen. Betrachten Sie zwei benachbarte Einträge als nach Suffix sortiertes Array.

      p     l
      v     v
   issippi
   ississippi
   

   p bezeichnet die Länge des längsten gemeinsamen Präfixes der beiden Einträge, l bezeichnet die Länge des zweiten Eintrags. Für ein solches Paar finden wir einen neuen unterschiedlichen Teilstring der Länge k für p < k ≤ l . Da p ≪ l häufig gilt, ist es unpraktisch, eine große Anzahl von Array-Einträgen für jedes Paar zu erhöhen. Stattdessen speichern wir das Array als Differenzarray, wobei jeder Eintrag k die Differenz zur Anzahl der k- mere zur Anzahl der ( k  - 1) -mers bezeichnet. Dies führt zu einer Aktualisierung des Formulars

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

   in eine viel schnellere Aktualisierung des Formulars

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

   Wenn wir sorgfältig beobachten, dass l immer unterschiedlich ist und tatsächlich jede 0 < l < n genau einmal erscheint, können wir die Subtraktionen weglassen und stattdessen 1 von jeder Differenz subtrahieren, wenn wir von Differenzen in Beträge umrechnen.

4. Differenzen werden in Beträge umgerechnet und ausgedruckt.

Der Quellcode

Diese Quelle verwendet libdivsufsort zum Sortieren von Suffix-Arrays. Der Code generiert beim Kompilieren mit diesem Aufruf eine Ausgabe gemäß der Spezifikation.

cc -O -o dsskmer dsskmer.c -ldivsufsort

Alternativ kann der Code eine Binärausgabe erzeugen, wenn er mit dem folgenden Aufruf kompiliert wird.

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

Um das höchste k zu begrenzen, für das k- mere gezählt werden sollen, geben Sie -DICAP = k an, wobei k die Grenze ist:

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

Kompilieren Sie mit, -O3wenn Ihr Compiler diese Option bereitstellt.

#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;
}

Das binäre Dateiformat kann mit folgendem Programm in das für Menschen lesbare Ausgabeformat konvertiert werden:

#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;
}

Beispielausgabe

Eine Beispielausgabe im Binärformat für die Datei chr22.fafinden Sie hier . Bitte bzip2 -dzuerst mit dekomprimieren . Die Ausgabe erfolgt nur im Binärformat, da sie viel besser komprimiert wird (3,5 kB gegenüber 260 MB) als die Ausgabe im lesbaren Format. Beachten Sie jedoch, dass die Referenzausgabe im unkomprimierten Zustand eine Größe von 924 MB hat. Möglicherweise möchten Sie eine Pipe wie die folgende verwenden:

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

C ++, k = 16, 37 Sekunden

Berechnet alle 16-mers in der Eingabe. Jedes 16-mer wird 4 Bit zu einem Symbol in ein 64-Bit-Wort gepackt (wobei ein Bitmuster für EOF reserviert ist). Die 64-Bit-Wörter werden dann sortiert. Die Antwort für jedes k kann abgelesen werden, indem man sich ansieht, wie oft sich die 4 * k oberen Bits der sortierten Wörter ändern.

Für den Testeingang verwende ich ca. 1,8 GB direkt unter dem Kabel.

Ich bin sicher, dass die Lesegeschwindigkeit verbessert werden könnte.

Ausgabe:

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

Programm:

#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);
    }
}

Kompilieren mit g++ -O3 kmer.cc -o kmerund ausführen mit ./kmer chr2.fa.

C ++ - eine Verbesserung gegenüber der FUZxxl-Lösung

Ich verdiene absolut keine Anerkennung für die Berechnungsmethode selbst, und wenn sich in der Zwischenzeit kein besserer Ansatz zeigt, sollte das Kopfgeld zu Recht an FUZxxl gehen.

#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;
}

Ich habe einfach Kasai et al. Algorithmus zur Berechnung von LCPs in O (n).
Der Rest ist eine bloße Anpassung des FUZxxl-Codes, wobei hier und da präzisere C ++ - Funktionen verwendet werden.

Ich habe den ursprünglichen Berechnungscode hinterlassen, um Vergleiche zu ermöglichen.

Da die langsamsten Prozesse die SA-Konstruktion und die LCP-Anzahl sind, habe ich die meisten anderen Optimierungen entfernt, um zu vermeiden, dass der Code für vernachlässigbare Gewinne überladen wird.

Ich habe die SA-Tabelle um das Präfix Null erweitert. Das erleichtert die LCP-Berechnung.

Ich habe keine Option zur Längenbeschränkung bereitgestellt. Der langsamste Prozess ist jetzt die SA-Berechnung, die nicht verkleinert werden kann (oder zumindest sehe ich nicht, wie es sein könnte).

Ich habe auch die Option für die binäre Ausgabe entfernt und die Anzeige auf die ersten 10 Werte beschränkt.
Ich gehe davon aus, dass dieser Code nur ein Proof of Concept ist, sodass keine überfüllten Displays oder gesättigten Festplatten erforderlich sind.

Erstellen der ausführbaren Datei

Ich musste das gesamte Projekt (einschließlich der Lite-Version vondivsufsort ) für x64 kompilieren , um das Win32 2Gb-Zuweisungslimit zu überwinden.

divsufsortCode gibt eine Reihe von Warnungen aus, da ints anstelle von size_ts häufig verwendet wird. Dies ist jedoch kein Problem für Eingaben unter 2 GB (für die ohnehin 26 GB RAM erforderlich wären: D).

Linux Build

Kompilieren main.cppund divsufsort.cVerwenden der Befehle:

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

Ich divsufsortgehe davon aus, dass die reguläre Bibliothek unter nativem Linux einwandfrei funktionieren sollte, solange Sie etwas mehr als 3 GB zuweisen können.

Aufführungen

Der Kasai-Algorithmus erfordert die inverse SA-Tabelle, die 4 weitere Bytes pro Zeichen für insgesamt 13 verbraucht (anstelle von 9 mit der FUZxxl-Methode).

Der Speicherverbrauch für den Referenzeingang liegt somit über 3 GB.

Andererseits wird die Rechenzeit dramatisch verbessert, und der gesamte Algorithmus befindet sich jetzt in 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.)

Weitere Verbesserungen

SA-Bau ist jetzt der langsamste Prozess.
Einige divsufsortTeile des Algorithmus sollen mit den mir unbekannten Funktionen eines Compilers parallelisiert werden. Falls erforderlich, sollte sich der Code jedoch leicht an klassischeres Multithreading anpassen lassen ( z. B. à la C ++ 11).
Die Bibliothek verfügt außerdem über eine Vielzahl von Parametern, darunter verschiedene Bucket-Größen und die Anzahl der unterschiedlichen Symbole in der Eingabezeichenfolge. Ich habe sie nur flüchtig betrachtet, aber ich vermute, dass das Komprimieren des Alphabets einen Versuch wert sein könnte, wenn Ihre Zeichenfolgen endlose ACTG-Permutationen sind ( und Sie verzweifelt nach Auftritten suchen ).

Es gibt auch einige parallelisierbare Methoden zur Berechnung von LCP aus SA, aber da der Code auf einem Prozessor, der etwas leistungsfähiger als mein mickriger [email protected] ist und der gesamte Algorithmus in O (n) ist, unter einer Minute laufen sollte , bezweifle ich dies wäre die Mühe wert.

C (kann auf meinem Computer bis zu 10 in einer Minute lösen)

Dies ist eine sehr einfache Lösung. Es konstruiert einen Baum der gefundenen k- mere und zählt sie. Um Speicherplatz zu sparen, werden Zeichen zuerst in Ganzzahlen von 0 bis n - 1 konvertiert, wobei n die Anzahl der verschiedenen Zeichen in der Eingabe ist. Daher müssen wir keinen Platz für Zeichen bereitstellen, die niemals angezeigt werden. Außerdem wird den Blättern weniger Speicher zugewiesen als anderen Knoten, um weiteren Speicher zu sparen. Diese Lösung benötigt zur Laufzeit auf meinem Computer etwa 200 MiB RAM. Ich verbessere es immer noch, vielleicht ist es in der Iteration sogar noch schneller!

Speichern Sie zum Kompilieren den folgenden Code in einer Datei mit dem Namen kmers.cund führen Sie ihn dann auf einem POSIX-ähnlichen Betriebssystem aus:

cc -O -o kmers kmers.c

Sie könnten ersetzen wollen -O3für , -Owenn Ihr Compiler unterstützt das. So führen Sie zunächst auspacken chr2.fa.gzin chr2.faund dann aus:

./kmers <chr2.fa

Dies erzeugt die Ausgabe Schritt für Schritt. Möglicherweise möchten Sie sowohl Zeit als auch Raum einschränken. Verwenden Sie so etwas wie

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

Ressourcen nach Bedarf zu reduzieren.

#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;
}

Verbesserungen

1. 8 → 9: Lesen Sie die gesamte Datei am Anfang, verarbeiten Sie sie einmal vor und behalten Sie sie im Kern. Dies verbessert den Durchsatz erheblich.
2. Verwenden Sie weniger Speicher, schreiben Sie die richtige Ausgabe.
3. Ausgabeformat erneut korrigieren.
4. Fix off-by-one.
5. 9 → 10: Werfen Sie nicht weg, was Sie bereits getan haben.