就地基数排序

显然,这是一种特定的DNA,而不是一般的,但它应该是快速的.

我很好奇这段代码是否真的有效,所以我在等待自己的生物信息学代码运行时测试/调试它.现在上面的版本实际上已经过测试和运行.对于每个5个碱基的1000万个序列,它比优化的内含子快约3倍.

I've never seen an in-place radix sort, and from the nature of the radix-sort I doubt that it is much faster than a out of place sort as long as the temporary array fits into memory.

Reason:

The sorting does a linear read on the input array, but all writes will be nearly random. From a certain N upwards this boils down to a cache miss per write. This cache miss is what slows down your algorithm. If it's in place or not will not change this effect.

I know that this will not answer your question directly, but if sorting is a bottleneck you may want to have a look at near sorting algorithms as a preprocessing step (the wiki-page on the soft-heap may get you started).

That could give a very nice cache locality boost. A text-book out-of-place radix sort will then perform better. The writes will still be nearly random but at least they will cluster around the same chunks of memory and as such increase the cache hit ratio.

I have no idea if it works out in practice though.

Btw: If you're dealing with DNA strings only: You can compress a char into two bits and pack your data quite a lot. This will cut down the memory requirement by factor four over a naiive representation. Addressing becomes more complex, but the ALU of your CPU has lots of time to spend during all the cache-misses anyway.

您当然可以通过以位为单位对序列进行编码来降低内存要求.您正在考虑排列,因此,对于长度为2,"ACGT"为16状态或4位.对于长度3,这是64个状态,可以以6位编码.因此,对于序列中的每个字母,它看起来像2位,或者像你说的那样,对于16个字符,大约32位.

如果有办法减少有效"单词"的数量,则可以进一步压缩.

因此,对于长度为3的序列,可以创建64个桶,可能是uint32或uint64.将它们初始化为零.迭代你非常大的3个字符序列列表,并按上述方式对它们进行编码.使用此作为下标,并增加该桶.
重复此操作,直到处理完所有序列.

接下来,重新生成列表.

按顺序迭代64个桶,对于在该桶中找到的计数,生成该桶所代表的序列的许多实例.
当所有桶都被迭代后,你就有了你的排序数组.

4的序列,加2位,所以会有256个桶.5的序列,加2位,所以会有1024个桶.

在某些时候,桶的数量将接近你的极限.如果从文件中读取序列,而不是将它们保留在内存中,则可以为存储桶提供更多内存.

我认为这比在原地进行排序要快,因为铲斗可能适合您的工作组.

这是一个显示该技术的黑客

#include 
#include 

#include 

using namespace std;

const int width = 3;
const int bucketCount = exp(width * log(4)) + 1;
      int *bucket = NULL;

const char charMap[4] = {'A', 'C', 'G', 'T'};

void setup
(
    void
)
{
    bucket = new int[bucketCount];
    memset(bucket, '\0', bucketCount * sizeof(bucket[0]));
}

void teardown
(
    void
)
{
    delete[] bucket;
}

void show
(
    int encoded
)
{
    int z;
    int y;
    int j;
    for (z = width - 1; z >= 0; z--)
    {
        int n = 1;
        for (y = 0; y < z; y++)
            n *= 4;

        j = encoded % n;
        encoded -= j;
        encoded /= n;
        cout << charMap[encoded];
        encoded = j;
    }

    cout << endl;
}

int main(void)
{
    // Sort this sequence
    const char *testSequence = "CAGCCCAAAGGGTTTAGACTTGGTGCGCAGCAGTTAAGATTGTTT";

    size_t testSequenceLength = strlen(testSequence);

    setup();


    // load the sequences into the buckets
    size_t z;
    for (z = 0; z < testSequenceLength; z += width)
    {
        int encoding = 0;

        size_t y;
        for (y = 0; y < width; y++)
        {
            encoding *= 4;

            switch (*(testSequence + z + y))
            {
                case 'A' : encoding += 0; break;
                case 'C' : encoding += 1; break;
                case 'G' : encoding += 2; break;
                case 'T' : encoding += 3; break;
                default  : abort();
            };
        }

        bucket[encoding]++;
    }

    /* show the sorted sequences */ 
    for (z = 0; z < bucketCount; z++)
    {
        while (bucket[z] > 0)
        {
            show(z);
            bucket[z]--;
        }
    }

    teardown();

    return 0;
}

sort(List elements, int prefix)
    if (elements.Count < THRESHOLD)
         return InMemoryRadixSort(elements, prefix)
    else
         return DiskBackedRadixSort(elements, prefix)

DiskBackedRadixSort(elements, prefix)
    DiskBackedBuffer[] buckets
    foreach (element in elements)
        buckets[element.MSB(prefix)].Add(element);

    List ret
    foreach (bucket in buckets)
        ret.Add(sort(bucket, prefix + 1))

    return ret

GATTACA