Menghitung k-mer

8

Tugasnya adalah untuk menghitung jumlah substring panjang k yang berbeda, untuk k = 1,2,3,4, .....

Keluaran

Anda harus output satu baris per kAnda berhasil menyelesaikan dengan satu nomor per baris output. Output Anda harus dalam urutan meningkat khingga Anda kehabisan waktu.

Skor

Skor Anda adalah k tertinggi yang dapat Anda peroleh di komputer saya dalam waktu kurang dari 1 menit.

Anda harus menggunakan http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/chr2.fa.gz sebagai masukan dan abaikan baris baru. Namun, kode Anda harus peka terhadap huruf besar-kecil.

Anda dapat mendekompresi input sebelum memulai timing.

Kode berikut (tidak efisien) menghitung jumlah 4-mers yang berbeda.

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

Batas memori

Untuk membuat kode Anda berperilaku baik di komputer saya dan untuk membuat tugas lebih menantang, saya akan membatasi RAM yang dapat Anda gunakan hingga 2GB dengan menggunakan

ulimit -v 2000000

sebelum menjalankan kode Anda. Saya sadar ini bukan cara solid untuk membatasi penggunaan RAM, jadi jangan gunakan cara imajinatif untuk mengatasi batasan ini dengan, misalnya, memunculkan proses baru. Tentu saja Anda dapat menulis kode multi-utas tetapi jika seseorang melakukannya, saya harus belajar cara membatasi total RAM yang digunakan untuk itu.

Tie Breaker

Dalam kasus dasi untuk beberapa maksimum ksaya akan menghitung berapa lama waktu yang dibutuhkan untuk memberikan hasil k+1dan yang tercepat menang. Jika mereka berjalan dalam waktu yang sama hingga dalam satu detik hingga k+1, pengiriman pertama menang.

Bahasa dan perpustakaan

Anda dapat menggunakan bahasa apa pun yang memiliki kompiler / juru bahasa / dll yang tersedia secara bebas. untuk Linux dan semua perpustakaan yang juga tersedia secara bebas untuk Linux.

Mesin saya Pengaturan waktu akan dijalankan di mesin saya. Ini adalah instalasi ubuntu standar pada Prosesor Delapan Core AMD FX-8350 pada Motherboard Asus M5A78L-M / USB3 (Socket AM3 +, 8GB DDR3). Ini juga berarti saya harus dapat menjalankan kode Anda. Sebagai konsekuensinya, hanya gunakan perangkat lunak gratis yang mudah tersedia dan harap sertakan instruksi lengkap cara menyusun dan menjalankan kode Anda.


Uji keluaran

Kode FUZxxl menampilkan yang berikut (tetapi tidak semua dalam 1 menit) yang saya yakini benar.

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

tabel liga

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

Seberapa banyak Anda dapat mengkhususkan kode Anda untuk input?

  • Jelas itu akan merusak kompetisi jika Anda hanya menghitung jawaban dan kode Anda menghasilkannya. Jangan lakukan itu.
  • Idealnya, apa pun yang perlu dipelajari kode Anda tentang data yang akan digunakan untuk berjalan lebih cepat, hal itu dapat dipelajari saat dijalankan.
  • Namun Anda dapat mengasumsikan input akan terlihat seperti data dalam file * .fa di http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/ .

sumber
Dalam J, solusi naif hanya akan menjadi `[: ~.]` Tetapi tebak itu tidak akan memotongnya.
FUZxxl
@ FuZxxl Yah ... seberapa besar aknya memberi Anda dalam satu menit dan banyak RAM yang digunakan?
Nah, untuk 2 ini kembali dalam waktu sekitar 5 detik, untuk 3 dibutuhkan 3,2 GB RAM dan membutuhkan waktu sekitar satu menit. Saya tidak mencoba apa yang akan dilakukan untuk empat. Biarkan saya mencoba ...
FUZxxl
@ FuZxxl Yah ... jangan ragu untuk mengirimkannya sebagai jawaban :) Ingat saya memotongnya pada 2GB RAM.
1
@Lembik saya tidak akan memberi tahu Anda. Saya saat ini sedang dalam proses mengubah solusi sehingga memberikan earlies output (tetapi secara keseluruhan lebih lambat). Butuh 9 menit untuk memproses data pada mesin saya.
FUZxxl

Jawaban:

7

C (≥ 4000)

Kode ini tidak menggunakan kurang dari 2 GB RAM (menggunakan sedikit lebih) atau menghasilkan setiap output dalam menit pertama. Tetapi jika Anda menunggu sekitar enam menit, ia mencetak semua jumlah k -mer sekaligus.

Opsi disertakan untuk membatasi k tertinggi yang kami hitung k -mers. Ketika k terbatas pada rentang yang masuk akal, kode berakhir dalam waktu kurang dari satu menit dan menggunakan kurang dari 2 GiB RAM. OP telah memberi peringkat solusi ini dan berakhir dalam waktu kurang dari satu menit untuk batas yang tidak lebih tinggi dari 4000.

Bagaimana cara kerjanya?

Algoritma ini memiliki empat langkah.

  1. Baca file input ke dalam buffer dan hapus baris baru.
  2. Akhiran-urutkan array indeks ke buffer input. Misalnya, sufiks string mississippiadalah:

    mississippi
    ississippi
    ssissippi
    sissippi
    issippi
    ssippi
    sippi
    ippi
    ppi
    pi
    i
    

    String ini diurutkan dalam urutan leksikografis adalah:

    i
    ippi
    issippi
    ississippi
    mississippi
    pi
    ppi
    sippi
    sissippi
    ssippi
    ssissippi
    

    Sangat mudah untuk melihat bahwa semua substring yang sama dengan panjang k untuk semua k ditemukan dalam entri yang berdekatan dari array yang diurutkan akhiran.

  3. Array integer diisi di mana kami menyimpan di setiap indeks k jumlah k -mers yang berbeda . Ini dilakukan dengan cara yang sedikit berbelit-belit untuk mempercepat proses. Pertimbangkan dua entri yang berdekatan dengan susunan suffix sort.

       p     l
       v     v
    issippi
    ississippi
    

    p menunjukkan panjang awalan umum terpanjang dari dua entri, l menunjukkan panjang entri kedua. Untuk pasangan seperti itu, kami menemukan substring panjang k yang berbeda untuk p < kl . Karena pl sering berlaku, tidak praktis untuk menambah sejumlah besar entri array untuk setiap pasangan. Sebagai gantinya kita menyimpan array sebagai array perbedaan, di mana setiap entri k menunjukkan perbedaan dengan jumlah k -mers ke jumlah ( k  - 1) -mers. Ini mengubah pembaruan formulir

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

    menjadi pembaruan formulir yang jauh lebih cepat

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

    Dengan mengamati dengan cermat bahwa l selalu berbeda dan pada kenyataannya setiap 0 < l < n akan muncul tepat satu kali, kita dapat menghilangkan pengurangan dan bukannya mengurangi 1 dari setiap perbedaan ketika mengkonversi dari perbedaan menjadi jumlah.

  4. Perbedaan dikonversi menjadi jumlah dan dicetak.

Kode sumber

Sumber ini menggunakan libdivsufsort untuk mengurutkan array sufiks. Kode menghasilkan output sesuai dengan spesifikasi saat dikompilasi dengan doa ini.

cc -O -o dsskmer dsskmer.c -ldivsufsort

alternatifnya kode dapat menghasilkan output biner ketika dikompilasi dengan doa berikut.

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

Untuk membatasi k tertinggi di mana k -mers dihitung, supply -DICAP = k di mana k adalah batasnya:

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

Kompilasi dengan -O3jika kompiler Anda menyediakan opsi ini.

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

Format file biner dapat dikonversi ke format output yang dapat dibaca manusia dengan program berikut:

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

output sampel

Contoh output dalam format biner untuk file chr22.fadapat ditemukan di sini . Silakan dekompresi bzip2 -dterlebih dahulu. Output disediakan dalam format biner hanya karena kompresnya jauh lebih baik (3,5 kB vs 260 MB) daripada output dalam format yang dapat dibaca manusia. Berhati-hatilah meskipun keluaran referensi memiliki ukuran 924 MB saat tidak dikompresi. Anda mungkin ingin menggunakan pipa seperti ini:

bzip2 -dc chr2.kmerdiffdat.bz2 | ./unbin | less
FUZxxl
sumber
Ini sangat menarik. Bisakah Anda merinci waktu yang diambil oleh masing-masing bagian? Tampaknya membangun susunan sufiks kurang dari 10% dari waktu.
1
langkah satu membutuhkan waktu sekitar tiga detik, langkah kedua membutuhkan waktu sekitar 20 detik, langkah ketiga mengambil seluruh istirahat. Waktu yang diambil oleh langkah empat dapat diabaikan.
FUZxxl
1
@FUZxxl apakah Anda memplot distribusi nilai p maksimum untuk melihat seberapa banyak Anda dapat mempercepat langkah 3 dengan memotong jika p> k (hanya menghitung untuk k-mer)? meskipun jika langkah 2 membutuhkan waktu 20 detik mungkin tidak ada banyak ruang di sana. penjelasan hebat btw
randomra
1
@Lembik Jangan gunakan cat. Gunakan pengalihan shell seperti di ./dsskmer <~/Downloads/chr2.fs. Kode perlu tahu berapa lama file input dan itu tidak mungkin dengan pipa.
FUZxxl
2
@Lembik Sepertinya tidak terlalu lambat untukku. Mungkin mereka hanya programmer yang buruk.
FUZxxl
7

C ++, k = 16, 37 detik

Menghitung semua 16-mers dalam input. Setiap 16-mer dikemas 4 bit ke simbol menjadi kata 64-bit (dengan pola satu bit dicadangkan untuk EOF). Kata-kata 64-bit kemudian disortir. Jawaban untuk setiap k dapat dibaca dengan melihat seberapa sering bit top 4 * k dari kata yang diurutkan berubah.

Untuk input tes saya menggunakan sekitar 1.8GB, tepat di bawah kawat.

Saya yakin kecepatan membaca dapat ditingkatkan.

Keluaran:

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

Program:

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

Kompilasi dengan g++ -O3 kmer.cc -o kmerdan jalankan ./kmer chr2.fa.

Keith Randall
sumber
Harap patuhi format keluaran; jangan k .
FUZxxl
Sebaliknya, ini adalah solusi keren.
FUZxxl
@ FuZxxl: diperbaiki.
Keith Randall
Berapa lama menurut Anda untuk menemukan panjang terpanjang dari setiap substring yang diulang? Saya pikir Anda secara otomatis mengetahui semua jawaban untuk k lebih lama dari itu.
@ Lembik: Substring berulang yang terpanjang dapat dilakukan dalam waktu linier. Tapi saya tidak berpikir itu sangat membantu - ada substring yang benar-benar panjang di input sampel, setidaknya ribuan simbol.
Keith Randall
4

C ++ - peningkatan dari solusi FUZxxl

Saya sama sekali tidak layak kredit untuk metode perhitungan itu sendiri, dan jika tidak ada pendekatan yang lebih baik muncul pada saat itu, karunia harus pergi ke FUZxxl dengan benar.

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

Saya hanya menggunakan Kasai et al. algoritma untuk menghitung LCP di O (n).
Sisanya hanyalah adaptasi dari kode FUZxxl, menggunakan fitur C ++ yang lebih ringkas di sana-sini.

Saya meninggalkan kode perhitungan asli untuk memungkinkan perbandingan.

Karena proses yang paling lambat adalah konstruksi SA dan jumlah LCP, saya menghapus sebagian besar optimasi lainnya untuk menghindari kekacauan kode untuk mendapatkan keuntungan yang dapat diabaikan.

Saya memperpanjang tabel SA untuk menyertakan awalan nol panjang. Itu membuat perhitungan LCP lebih mudah.

Saya tidak memberikan opsi batasan panjang, proses paling lambat saat ini adalah perhitungan SA yang tidak dapat dirampingkan (atau setidaknya saya tidak melihat bagaimana itu bisa dilakukan).

Saya juga menghapus opsi output biner dan tampilan terbatas pada 10 nilai pertama.
Saya menganggap kode ini hanya bukti konsep, jadi tidak perlu mengacaukan tampilan atau jenuh disk.

Membangun executable

Saya harus mengkompilasi seluruh proyek (termasuk versi litedivsufsort ) untuk x64 untuk mengatasi batas alokasi Win32 2Gb.

divsufsortkode melempar banyak peringatan karena penggunaan ints, bukan size_ts, tetapi itu tidak akan menjadi masalah untuk input di bawah 2Gb (yang tetap membutuhkan 26Gb RAM: D).

Linux membangun

kompilasi main.cppdan divsufsort.cgunakan perintah:

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

Saya menganggap divsufsortperpustakaan biasa harus bekerja dengan baik di Linux asli, selama Anda dapat mengalokasikan sedikit lebih dari 3Gb.

Pertunjukan

Algoritma Kasai membutuhkan tabel SA terbalik, yang memakan hingga 4 byte lebih per karakter untuk total 13 (bukan 9 dengan metode FUZxxl).

Konsumsi memori untuk input referensi dengan demikian di atas 3Gb.

Di sisi lain, waktu perhitungan ditingkatkan secara dramatis, dan keseluruhan algoritma sekarang dalam 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.)

Perbaikan lebih lanjut

Konstruksi SA sekarang merupakan proses paling lambat.
Beberapa bit dari divsufsortalgoritma dimaksudkan untuk diparalelkan dengan fitur bawaan apa pun dari kompiler yang tidak saya ketahui, tetapi jika perlu kodenya harus mudah beradaptasi dengan multithreading yang lebih klasik ( à la C ++ 11, misalnya).
Lib juga memiliki satu truk penuh parameter, termasuk berbagai ukuran ember dan jumlah simbol yang berbeda dalam string input. Saya hanya melihat sekilas pada mereka, tapi saya curiga mengompresi alfabet mungkin patut dicoba jika string Anda adalah permutasi ACTG yang tak ada habisnya ( dan Anda sangat membutuhkan pertunjukan).

Ada beberapa metode yang dapat diparalelkan untuk menghitung LCP dari SA juga, tetapi karena kode harus berjalan di bawah satu menit pada prosesor yang sedikit lebih kuat daripada [email protected] saya yang kecil dan seluruh algoritme dalam O (n), saya ragu ini akan sepadan dengan usaha.


sumber
2
+1 untuk memberi kredit saat kredit jatuh tempo;) (dan percepatan bagus!)
Martin Ender
1
Ini rapi. Perbaikan yang bagus.
FUZxxl
Saya bahkan tidak tahu tentang array LCP. Kasai et al. Algoritma juga sangat rapi. Agak dikatakan bahwa itu membutuhkan banyak memori.
FUZxxl
@ FuZxxl ah well, itu selalu sama tradeoff kecepatan / memori yang sama, tapi saya pikir mem 45%. kenaikan biaya adalah harga yang dapat diterima untuk mencapai kompleksitas O (n).
2
@ kuroineko Saya punya ide bagaimana membangun data yang kita butuhkan dalam waktu linier tanpa menggunakan array LCP. Biarkan saya periksa ...
FUZxxl
1

C (dapat menyelesaikan hingga 10 dalam satu menit di mesin saya)

Ini adalah solusi yang sangat sederhana. Itu membangun pohon k -mers yang ditemukan dan menghitungnya. Untuk menghemat memori, karakter terlebih dahulu dikonversi menjadi bilangan bulat dari 0 ke n - 1 dengan n adalah jumlah karakter yang berbeda dalam input, jadi kami tidak perlu menyediakan ruang untuk karakter yang tidak pernah muncul. Selain itu, lebih sedikit memori yang dialokasikan untuk daun daripada untuk node lain untuk menghemat memori lebih lanjut. Solusi ini menggunakan sekitar 200 MiB RAM saat runtime di komputer saya. Saya masih memperbaikinya, jadi mungkin dalam iterasi mungkin lebih cepat!

Untuk mengkompilasi, simpan kode di bawah ini dalam file bernama kmers.cdan kemudian jalankan pada sistem operasi seperti POSIX:

cc -O -o kmers kmers.c

Anda mungkin ingin mengganti -O3untuk -Ojika dukungan compiler Anda yang. Untuk menjalankan, pertama buka paket chr2.fa.gzke chr2.fadan kemudian jalankan:

./kmers <chr2.fa

Ini menghasilkan langkah demi langkah keluaran. Anda mungkin ingin membatasi waktu dan ruang. Gunakan sesuatu seperti

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

untuk mengurangi sumber daya sesuai kebutuhan.

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

Perbaikan

  1. 8 → 9: Baca seluruh file di awal, pra-proses sekali dan simpan di-core. Ini sangat meningkatkan throughput.
  2. Gunakan lebih sedikit memori, tulis output yang benar.
  3. Perbaiki format output lagi.
  4. Memperbaiki satu per satu.
  5. 9 → 10: Jangan membuang apa yang sudah kamu lakukan.
FUZxxl
sumber
Outputnya harus benar-benar satu angka per baris. Tampaknya saat ini untuk mengeluarkan string dari file input.
@Lembik Ah, saya mengerti! Sepertinya saya salah mengerti format output. Beri aku waktu untuk memperbaikinya.
FUZxxl
Format Output @Lembik diperbaiki lagi. Jika Anda suka, saya dapat menambahkan kode yang membunuh program setelah satu menit.
FUZxxl
Terima kasih. Anda saat ini adalah pemenang :) timeout 60sberfungsi OK untuk saya sehingga tidak perlu membangun cara untuk membunuh kode setelah 1 menit.
Anda juga bisa menggunakannya ulimit -t 60.
FUZxxl