65.9K
CodeProject 正在变化。 阅读更多。
Home

使用C++在视频内存中进行FASTA核苷酸序列缓存

tugrulGtx
starIconstarIconstarIconstarIconstarIcon

5.00/5 (6投票s)

2021 年 2 月 25 日

GPL3

15分钟阅读

viewsIcon

6221

通过索引访问缓存于显存中的 FASTA 格式文件（*.fna, *.faa）的核苷酸序列

引言

当一个应用程序需要构建一个可随机访问的碱基对序列列表，而这需要一个庞大的核苷酸序列文件（例如 FASTA 格式）时，它可能需要比系统可用内存（RAM）更多的空间。有时，系统没有空的内存插槽来升级，以便在大型数据集上测试核苷酸序列匹配算法原型。但是，有些库可以直接从磁盘驱动器索引相关数据，而无需消耗可用的内存。另一方面，这个工具借助霍夫曼编码算法和 OpenCL 数据传输，在显卡内存上做了类似的事情。它非常适用于 RAM 较低 + VRAM 较高，或者廉价 HDD + 廉价 GPU + 昂贵 CPU 的非对称系统。

背景

该序列缓存算法的主要部分包括：

  • 对序列和描述符进行完全数据压缩
  • 将数据分布到多个显卡以增加带宽和容量
  • 在访问序列或描述符的索引时进行部分解压

压缩算法是霍夫曼编码。首先，解析所有数据以统计每个符号出现的次数。对于生物信息学来说，这至少包括 'A'、'C'、'G' 和 'T'。有些文件有更多的符号，但仅仅 4 个字母就能在压缩后使其大小减少到约 1/4。

统计完成后（这是第一遍扫描），会构建一棵霍夫曼树。构建过程如下：

  • 按出现次数的升序对所有计数进行排序
  • 取前两个出现次数，并将它们放入一个新的出现“节点”（该节点还有一些其他成员，用于在编码和解码时进行簿记）
  • 将这个新节点放回同一个数组中（使用优先队列代替数组也可以）
  • 再次对数组进行排序
  • 重复此过程，直到数组中只剩下一个节点

在 std::vector 上构建霍夫曼树的实际代码

        while(sortedNodes.size()>1)
        {
            // get "least count" nodes
            Node<T> node1 = sortedNodes[0];
            Node<T> node2 = sortedNodes[1];

            // prepare a parent node for the two nodes
            Node<T> newNode;
            newNode.count = node1.count + node2.count;
            newNode.data=0;
            newNode.leaf1 = node1.self;
            newNode.leaf2 = node2.self;
            newNode.self = ctr;
            newNode.isLeaf=false;

            // remove the two nodes
            sortedNodes.erase(sortedNodes.begin());
            sortedNodes.erase(sortedNodes.begin());

            // add new node
            sortedNodes.push_back(newNode);

            // also add for referencing from decoding method
            referenceVec.push_back(newNode);

            // sort the vector so that "least count" nodes go to lower indices
            std::sort(sortedNodes.begin(), sortedNodes.end(),
            [](const Node<T> & n1, const Node<T> & n2){ return n1.count<n2.count;});
            ctr++;
        }

        root = sortedNodes[0];

这个结果节点元素被称为根节点，并被用作解码某个符号的路径入口点。

template<typename T>
struct Node
{
    size_t count;
    int self;
    int leaf1;
    int leaf2;
    T data;
    bool isLeaf;
};

这使得从二进制编码的数据路径中找到一个符号变得容易。根节点没有位。但 leaf1 被命名为 '0'，leaf2 被命名为 '1'。因此，当读取压缩数据时，如果第一个元素是 '1'，则选择 leaf2 节点索引。之后，从压缩数据中读取下一个位。然后，再次选择 leaf1 或 leaf2 路径，直到 isLeaf 为 true。当 isLeaf 为 true 时，所选节点是路径中的最后一个节点，并在其数据成员中包含解码后的符号（如 'A', 'C', ..）。

每个符号的“路径”由一个递归函数计算，该函数在树遍历的每个深度都持有一个临时的 std::vector<bool> 变量。每个新的深度都会向向量中添加一个新的位，并将结果向量赋给另一个数据结构，以便在解码期间访问。

通常，仅凭 count、leaf1 和 leaf2 就足以同时进行编码和解码。但为了优化速度，isLeaf 成员被用作快捷方式，而不是在解码时同时检查 leaf1 和 leaf2。数据类型是 unsigned char，这意味着这里进行的是字节级压缩。为了提高压缩率，未来也许可以使用字节包，如 short、int 甚至是 size_t。但提高压缩分辨率也会增加不适合压缩的数据的树大小。例如，如果一个 FASTA 文件包含的字母不止“ACGT”，那么对于基于 int 的压缩，霍夫曼树中的叶节点数量可能增长到 256 * 256 * 256 * 256 个。这对于缓存访问不利，也违背了压缩的目的（通过计算工作来节省带宽和内存容量）。目前，仅使用 unsigned char，霍夫曼树最多包含几千字节的节点数据。

在文件的下一遍扫描中（第二遍），会推导出其总比特数。由于内存可能不足以容纳大文件，需要额外一遍扫描来计算压缩数据的比特数。FASTA 文件的描述符和序列都被压缩，因此有一个大的比特流同时包含序列比特流和描述符比特流。但两者都有各自独特的霍夫曼树，这样一来，虽然描述符通常有 256 个字符且压缩效果较差，但序列部分仍然可以根据其自身的符号表得到很好的压缩（例如，仅有 A 和 B 两个字母可以被压缩为 1 和 0，实现 8 倍压缩）。为了计算比特数，每个字母都映射到一个比特流大小。例如，如果腺嘌呤核碱基（'A' 符号）的霍夫曼树路径是 '1 - 1 - 0'，那么它的比特流就有 3 个元素。'A' 的比特数就增加 3。

一旦知道比特数，就会分配一个由显存支持的虚拟数组（该数组具有静态大小和对显卡的交错访问，作为某种“分页”机制，可对任何数据进行线程安全访问，并利用多线程访问中 PCIe 的高带宽优势）。（其工作原理：VirtualMultiArray - C++ 虚拟数组实现）

如果虚拟数组分配成功，则开始第三遍扫描。这一次，不再是计算比特数，而是将这些比特复制到虚拟数组中。每个比特都被放入一个 unsigned char 元素的相关比特位置。第一个比特进入第一个字节的第一个比特位，第二个比特进入第一个字节的第二个比特位……直到最后一个比特进入最后一个字节的第 k 个比特位，此过程完成。

读取一个大文件 3 次速度很慢。这仅用于初始化。一旦完成，访问一个序列只需要微秒级的 PCIe 数据复制和微秒级的霍夫曼解码（并且时间随所需序列的大小线性增加）。

每次访问一个序列索引都会经过显存缓存的“后备存储”、内存缓存的“页面”和解码（小序列在 L1 缓存中，大序列在内存中）。对于单线程访问，这些是顺序任务。在多线程访问中，这些延迟相互隐藏，因此当一个线程从显存加载数据时，另一个线程可以从内存加载数据，而另一个线程则在进行霍夫曼解码计算。

Using the Code

要启动一个从 FASTA 格式文件中读取序列的最小示例，从仓库（https://github.com/tugrul512bit/FastaGeneIndexer）获取头文件后，至少需要准备几件事。

首先，内部的虚拟数组实现（VirtualMultiArray, https://github.com/tugrul512bit/VirtualMultiArray）使用 OpenCL。这需要链接器链接到其驱动程序库（对于 Ubuntu+Nvidia 是 /usr/local/cuda/lib64 文件夹，Windows 是 system32 文件夹），并包含相关的头文件目录以使用 OpenCL 数据复制功能。

其次，它在其中一个解析阶段使用了 OpenMP。这需要为编译器和链接器启用相关的 OpenMP 标志。一些编译器已经内置了此功能。

最后，它使用 C++17 获取文件大小，其他部分使用 C++14，并且项目需要是 64 位（在 g++10 中除非另有说明，否则这是默认设置，也适用于 g++7.5 和 msvc 2019）。

在开发计算机（Ubuntu 18.04, Nvidia Nsight, g++10）上，编译和链接命令如下：

g++-10 -std=c++17 -O3 -g3 -Wall -c -fmessage-length=0 -mavx -march=native 
-funroll-loops -mtune=native -funsafe-math-optimizations -ftree-vectorize 
-ffast-math  -fopenmp -pthread -fPIC -MMD -MP -MF"src/main.d" -MT"src/main.d" 
-o "src/main.o" "../src/main.cpp"
g++-10 -L/usr/local/cuda/lib64 -o "Virtualizer2"  ./src/main.o   -lOpenCL -lgomp -lpthread

除了需要 VirtualMultiArray 头文件外，其他所有功能都包含在单个头文件中。

#include"FastaGeneIndexer.h"

读取一个 FASTA 文件只需一行代码。

FastaGeneIndexer cache("./data/influenza.fna");

如果文件解析通过并且需要查看霍夫曼树的内容，可以提供一个调试标志。

bool debug = true;
FastaGeneIndexer cache("./data/influenza.fna", debug);

如果虚拟数组/OpenCL API 的参数无效，它会抛出异常，可以像这样捕获它们：

try
{
   bool debug = true;
   FastaGeneIndexer cache("./data/influenza.fna", debug);

}
catch(std::exception & e)
{
   std::cout<< e.what() <<std::endl;
}

然后就可以读取序列了。

#include "FastaGeneIndexer.h"
int main(int argC, char** argV)
{
    try
    {
        bool debug = true;
        FastaGeneIndexer cache("./data/influenza.fna", debug);
        std::string sequence = cache.getSequence(5);
        std::string descriptor = cache.getDescriptor(5);

        std::cout<< sequence <<std::endl;
        std::cout<< descriptor <<std::endl;
    }
    catch(std::exception & e)
    {
        std::cout<< e.what() <<std::endl;
    }
    return 0;
}
输出
ATGAGTCTTCTAACCGAGGTCGAAACGTACGTTCTCTCTATCGTCCCGTCAGGCCCCCTCAAAGCCGAGATCGCACAGAGACTT
GAAGATGTCTTTGCAGGGAAGAACACCGATCTTGAGGCTCTCATGGAATGGCTAAAGACAAGACCAATCCTGTCACCTCTGACT
AAGGGGATTTTAGGATTTGTGTTCACGCTCACCGTGCCCAGTGAGCGGGGACTGCAGCGTAGACGCTTTGTCCAAAATGCTCTT
AATGGGAACGGAGATCCAAATAACATGGACAAAGCAGTTAAACTGTATAGGAAGCTTAAGAGGGAGATAACATTCCATGGGGCC
AAAGAAATAGCACTCAGTTATTCTGCTGGTGCACTTGCCAGTTGTATGGGCCTCATATACAACAGGATGGGGGCTGTGACCACT
GAAGTGGCATTTGGCCTGGTATGCGCAACCTGTGAACAGATTGCTGACTCCCAGCATCGGTCTCATAGGCAAATGGTGACAACA
ACCAATCCACTAATCAGACATGAGAACAGAATGGTTCTAGCCAGCACTACAGCTAAGGCTATGGAGCAAATGGCTGGATCGAGT
GAGCAAGCAGCAGAGGCCATGGATATTGCTAGTCAGGCCAGGCAAATGGTGCAGGCGATGAGAACCATTGGGACTCATCCTAGC
TCCAGTGCTGGTCTAAAAGATGATCTTCTTGAAAATTTGCAGGCCTATCAGAAACGAATGGGGGTGCAGATGCAACGATTCAAG
TGATCCTCTCGTTATTGCAGCAAATATCATTGGGATCTTGCACTTGATATTGTGGATTCTTGATCGTCTTTTTTTCAAATGCAT
TTATCGTCGCTTTAAATACGGTTTGAAAAGAGGGCCTTCTACGGAAGGAGTGCCAGAGTCTATGAGGGAAGAATATCGAAAGGA
ACAGCAGAATGCTGTGGATGTTGACGATGGTCATTTTGTCAACATAGAGCTGGAGTAA
gi|60460|gb|X08089|Influenza A virus (A/WS/7/1933(H1N1)) 
ORF1 for M1 protein and ORF2 for M2 protein, genomic RNA

FastaGeneIndexer 构造函数接受一个文件名字符串参数，并初始化访问描述符和序列所需的一切。解析文件、构建霍夫曼树、分配虚拟数组，所有工作都在构造函数中完成。对于大型输入文件，进行多次扫描需要一些时间，因此启用调试有助于查看压缩准备工作是否有进展。

无需显式释放任何资源。内部数据结构的析构函数会自动释放所有资源，如 OpenCL 缓冲区、页面和霍夫曼节点。

最大带宽在很大程度上依赖于多线程访问、多个显卡的存在以及处理器的单线程性能。

对单线程首次访问的简单基准测试结果：

#include "FastaGeneIndexer.h"
#include "CpuBenchmarker.h"
int main(int argC, char** argV)
{
    try
    {
        bool debug = true;
        FastaGeneIndexer cache("./data/influenza.fna", debug);

        {
            CpuBenchmarker bench(0,"read 1");
            cache.getSequence(10);
        }
    }
    catch(std::exception & e)
    {
        std::cout<< e.what() <<std::endl;
    }
    return 0;
}
输出
read 1: 104912 nanoseconds

这主要是因为首次执行 OpenCL 数据复制（导致驱动程序在某处惰性初始化一些内部状态），以及将数据从显存读入页面（即内存缓存）。下一次读取性能会更好。

        {
            CpuBenchmarker bench(0,"read 1");
            cache.getSequence(10);
        }
        
        {
            CpuBenchmarker bench(0,"read 2 same location");
            cache.getSequence(10);
        }
输出
read 1: 148684 nanoseconds    
read 2 same location: 77568 nanoseconds

仅两次访问可能无法展示真实世界的性能。CPU 也需要处于较高频率状态才能获得更好的基准测试结果。

        // heating cpu
        for(int i=0;i<10000;i++)
        {
            cache.getSequence(i);
        }

        {
            CpuBenchmarker bench(0,"read 10k sequences, single thread",10000);
            for(int i=0;i<10000;i++)
            {
                cache.getSequence(i);
            }
        }
输出
read 1: 145723 nanoseconds    
read 2 same location: 78204 nanoseconds    
read 10k sequences, single thread: 300136627 nanoseconds 
(throughput = 30013.66 nanoseconds per iteration)

这是在 FX8150@2.1GHz + 单通道 DDR3 1333MHz 内存 + (PCIe v2.0 x8 + x4 + x4) 上的结果。CPU 的核心性能和内存延迟越好，解码延迟就越低。解码意味着从显卡获取数据，在霍夫曼树上进行一些树遍历，以及在容器之间多次移动数据。

主要思想是对缓存进行并发访问，以获得比老式 HDD 更高的带宽。

       // heating cpu
        for(int i=0;i<10000;i++)
        {
            cache.getSequence(i);
        }

        size_t count=0;
        {
            CpuBenchmarker bench(0,"read 10k sequences, single thread",10000);
            for(int i=0;i<10000;i++)
            {
                count += cache.getSequence(i).size();
            }
        }

        {
            CpuBenchmarker bench(count,"read 10k sequences, multi-thread",10000);

            #pragma omp parallel for
            for(int i=0;i<10000;i++)
            {
                cache.getSequence(i);
            }
        }
输出
read 1: 146718 nanoseconds    
read 2 same location: 78615 nanoseconds    
read 10k sequences, single thread: 302150806 nanoseconds  
(throughput = 30215.08 nanoseconds per iteration)

read 10k sequences, multi-thread: 55878726 nanoseconds 
(bandwidth = 203.29 MB/s)      (throughput = 5587.87 nanoseconds per iteration)

每次读取的平均时间扩展性低于核心数（8），但大于核心模块数（4），并且单通道 1333MHz 内存对此有一定影响。为了更清楚地看到这一点，将 CPU 频率设置回默认值（3.6GHz）（经过 47 秒的初始化后）：

read 1: 154930 nanoseconds    
read 2 same location: 79527 nanoseconds    
read 10k sequences, single thread: 177180777 nanoseconds     
(throughput = 17718.08 nanoseconds per iteration)

read 10k sequences, multi-thread: 38953637 nanoseconds     
(bandwidth = 291.62 MB/s)      (throughput = 3895.36 nanoseconds per iteration)

CPU 频率提高 71%，平均访问时间提高了 43%。单线程访问延迟提高了 70%。这非常接近由单通道 DDR3 1333MHz 内存引起的内存带宽瓶颈。内存带宽由 PCI-E 数据传输和解码共享。解码序列过长意味着在 getSequence() 中容器之间移动的数据更多，留给 PCI-E 传输的带宽就更少。因此，对于某些只有 5-10 个序列（但每个序列都长达数兆字节）的 FASTA 文件，多线程的优势不大（对于当前版本），但性能仍优于单线程。

不同的系统会有不同的瓶颈。Threadripper 可能在内存上遇到瓶颈，而 Epyc 在使用 1-2 个显卡时可能会遇到 PCIe 瓶颈。

测试期间，示例文件为 1.4GB，仅使用了系统显存的 15%。

特点

子序列查询

用户也可以查询一个序列中 K 个核苷酸。

这个重载方法（FastaGeneIndexer::getSequence(index) 的重载）接受两个额外的参数：子序列在所选序列中的起始位置和其长度。

FastaGeneIndexer cache("./data/yersinia_pestis_genomic.fna");

// getting 100 nucleotide symbols from fourth(index=3) sequence,
// starting at its 101th (index=100) symbol
std::cout << cache.getSequence(3,100,100) << std::endl;
输出
GACAGTTATGGAAATTAAAATCCTGCACAAGCAGGGAATGAGTAGCCGGGCGATTGCCAGAGAACTGGGGAT
CTCCCGCAATACCGTTAAACGTTATTTG

它不会带有换行符地写入。所有存储在显卡中的数据都是纯粹的核苷酸名称和描述符字符串。

由于此方法不获取完整的序列数据，它完成得更快。在 3.6GHz 的 FX8150 上单符号查询的性能：

sub-sequence query (sequential-access latency benchmark, single thread):
5603564809 nanoseconds     
(throughput = 1160.19 nanoseconds per iteration)

sub-sequence query (random-access latency benchmark, single thread):
5432723756 nanoseconds     
(throughput = 1167.39 nanoseconds per iteration)

使用多线程进行 1000 符号查询的性能：

bandwidth benchmark, sequential access, multi-threaded:
166062 nanoseconds
(bandwidth = 423.37 MB/s) (throughput = 2372.31 nanoseconds per iteration)

(just 150 microseconds of benchmarking contains some noise)
bandwidth benchmark, random access, multi-threaded:
157730 nanoseconds
(bandwidth = 445.73 MB/s) (throughput = 2253.29 nanoseconds per iteration)

通过描述符名称查询序列和子序列

调用一次

FastaGeneIndexer::initDescriptorIndexMapping()

方法后，用户可以通过其描述符名称获取（子/完整）序列。

此方法接受一个字符串参数作为描述符名称，其余参数与子序列查询相同。

<code>FastaGeneIndexer::getSequenceByDescriptor(std::string, size_t, size_t)</code>
bool debug = true;
FastaGeneIndexer cache(<span class="pl-pds">"./data/influenza.fna"</span>, debug);
cache.initDescriptorIndexMapping();

// use descriptor string to query a sub-sequence of 50 length 
// at 50th(0-based)position
std::cout << cache.getSequenceByDescriptor(cache.getDescriptor(1),50,50) << std::endl;
输出
CGTACGTTCTCTCTATCGTCCCGTCAGGCCCCCTCAAAGCCGAGATCGCA

整个文件的总符号数

用户可以获取预先计算的符号计数值。

FastaGeneIndexer::getSymbolCount(unsigned char)
bool debug = true;
FastaGeneIndexer cache("./data/influenza.fna", debug);
std::cout << cache.getSymbolCount('A') << std::endl;
std::cout << cache.getSymbolCount('C') << std::endl;
std::cout << cache.getSymbolCount('G') << std::endl;
std::cout << cache.getSymbolCount('T') << std::endl;
输出
437637159
249039370
309033962
304296666

已索引的序列总数

方法 FastaGeneIndexer::n() 返回已索引（并已缓存）的序列数量。例如，influenza.fna 这个 1.4GB 的文件中有超过 80 万个序列。

并行序列查询

FastaGeneIndexer::getSequenceParallel(size_t index) 专为快速获取非常长的序列（例如人类染色体每个序列超过 250MB）而优化。它在底层使用 OpenMP 执行一个非常简单的 for 循环，并行获取数据块并将它们连接成一个 string 结果。由于解码需要时间，它利用所有核心来获得良好的加速效果。将工作分解成更小的块（大小为 pageSize，对于大文件是 256k）可以利用缓存并提高整体读取性能。

可选的额外内存使用

将 FastaGeneIndexer 类的第三个构造函数参数设置为 true 时，会启用使用更多内存的缓存，以允许多核 CPU 能够重叠序列查询，从而提高整体吞吐量。

可选择限制第一张显卡的显存使用

在多 GPU 系统中，第一张/主显卡总是被操作系统使用，会损失约 400MB（在 2GB 显卡上）的显存用于操作系统功能。因此，为 FastaGeneIndexer 构造函数添加了第四个参数。当此参数为 true 时，缓存会在主显卡上分配少 25% 的显存空间，并让其他显卡完全填满其显存与霍夫曼编码数据。例如，当在一个 3 x 2GB 显卡的系统上缓存 4 个 3.2GB 的 chromosome.fa 文件时，启用此参数会填充第一张显卡的 85%，而不是溢出。另外两张显卡则可毫无问题地使用其 95% 的显存。

如何查找多个 FASTA 文件之间的共同序列？

这就像在两个整数向量之间寻找交集，只是需要额外记录每个序列的长度。在这里，我们寻找完全相同的核碱基序列（每个序列都有自己的描述符，但本例不考虑描述符数据，只考虑序列部分），涉及 4 个文件，每个文件 3.2GB，包含 194 个序列（总共 194 x 4 个序列，总计 12.8GB）。

算法

  • 对每个文件，独立地按其长度对所有序列索引进行排序。
  • 选择前两个已排序的索引数组。
  • 执行“交集算法”
    • while ( 两个数组的计数器都在范围内 )
    • {
      • 找到一个两个序列长度相等的点
        • 扫描两个数组中整个相同长度的块
        • 暴力比较（对两个块的所有元素）
          • 从数组 1 的块中获取 x
          • 从数组 2 的块中获取 y
          • compare
            • 如果相等，则写入结果 map (std::map)
    • }
  • 重复 N-1 次，但这次使用包含交集元素的“结果 (std::map)”和下一个文件进行比较。
    • 也使用 result2 map 作为输出，以保留“result”的旧值
    • 在 result 和 result2 之间进行乒乓操作以保留中间结果
  • 当所有文件都比较完毕后，根据文件数量是偶数还是奇数，查看 result 或 result2 数组。

您可以在这里找到测试函数代码：https://github.com/tugrul512bit/FastaGeneIndexer/blob/main/test_find_common_sequences.cpp

从链接中的代码可以看出，排序实际上不是原地排序。这个工具只是一个只读缓存。排序是在数组之外通过对索引进行排序完成的，没有移动任何一个数据字节。只是在第二层间接寻址（std::vector, std::map）上移动索引。在这里，std::map 对于处理数组之间的哈希冲突非常有用，因为它无需任何检查即可自动删除重复项。

对于只有 4GB 内存、6GB（组合）显存和 FX8150（3.6GHz） 的开发计算机，基准测试结果如下：

  • 初始化时间 = 181 秒 = 70 MB/s 平均 HDD 读取 + 霍夫曼编码 吞吐量
    • 仅解析 4 个文件，导致只能使用 4 个线程
    • 部分受限于老旧 HDD 最高约 150 MB/s 的读取吞吐量（仅限顺序读取）
    • 部分由 Ubuntu 自身的文件缓存加速
    • 主要是编码瓶颈
  • 计算时间 = 72 秒 = 355 MB/s 平均显存数据获取 + 霍夫曼解码 + 比较吞吐量
    • 使用 CPU 的所有线程
    • 主要是解码瓶颈
    • 显卡越多越好
    • 字符串比较仅使用 ==，它由编译器优化，完成得很快
    • 假设 355 MB/s 的值是基于线性扫描。如果发生了大量暴力比较，那么实际带宽会更高
    • 系统内存是单通道 1333MHz = 10GB/s = 这是基于内存扫描且几乎不做任何有用工作的绝对最大带宽。一旦 std::string 数据在容器之间复制，带宽就会损失。查询方法中的 const 正确性有助于优化这部分。
  • CPU：高负载（计算时使用率 85%+）
  • 系统：完全响应
    • 这就是主要思想

关注点

将霍夫曼树作为节点和指针的组合来维护，比使用节点向量和整数索引要慢。这可能是因为指针的大小（系统中为 8 字节）是整数索引的两倍，并且会跳转到远离缓存的位置。即使是 16 位索引也可以使用，因为这只是字符级压缩，任何时候都不需要超过 64k 个节点。但是在压缩流中逐位跟踪时，32 位整数（或任何适合系统的大小）就足够好了。Leaf1 - leaf2 索引指向节点向量的元素。这有助于缓存访问节点，而无需跳转到内存中的遥远位置。

经过速度相关的优化后，std::map（用于统计符号出现次数）成为前两次文件解析过程中的瓶颈。与其像下面这样映射：

std::map<unsigned char, size_t> counts;
if(counts.find(symbol)==counts.end())
{
    counts[symbol]=1;
}
else
{
    counts[symbol]++;
}

使用简单的数组要快得多。

unsigned char counts[256];

counts[symbol]++

同样，将每个符号的压缩位（std::vector<bool>）保存在 std::map 中，比使用一个大小为 256 的简单数组要慢得多。稍微增加内存消耗带来了更好的解析性能。

学习霍夫曼编码是一次很好的经历，对于文本压缩，即使是基于字节的，它也非常出色。对于仅含“ACTG”的组合，它能实现高达约 3.9 倍的压缩。对于包含完整字母表和一些其他符号的情况，它仍然能将约 600MB 压缩到约 450MB，这比什么都不做要好得多。

如果 FASTA 格式文件的“序列”部分包含单词（而不是字母），那么通过对算法进行单词级调整，可以获得更好的压缩效果。可能会出现类似 { 0: "hello", 1:"world" } 的情况，实现 40 倍的压缩。但只有描述符部分有单词，而且描述符比序列小得多，所以没有针对它们进行基于单词的霍夫曼编码优化，而是简单地将它们（但使用它们自己的霍夫曼树）连接到与序列相同的比特流中。

在当前版本中，主要瓶颈是：

  • 初始化（3 次文件解析，最后一次由于霍夫曼编码而耗时更长）
  • CPU 单线程性能（用于使用霍夫曼树进行解码）
  • 内存带宽（由显存数据传输和解码共享）
  • PCI-E 带宽（仅当其他部分完全优化或有足够的硬件规格支持时）

由于这些原因，它仅适用于长时间对序列进行大量多线程访问。也许将来可以优化以采用更好的压缩方案，例如查找最可能出现的核苷酸重复的游程长度（如 AAACCCCAAAGTAAACCCCAAAGT），并将它们视为单个符号。一旦压缩率提高，PCIe 或内存带宽就变得不那么重要，CPU 核心性能也能得到更好的利用（前提是解码后的数据能放入缓存）。重复的序列也可以被视为一个符号。但查找重复序列会在初始化中增加另一层处理。

历史

  • 2021年2月24日：添加了基本信息
  • 2021年2月25日：改进了文章展示，加入了所用算法的说明
  • 2021年2月27日：添加了新功能。子序列查询、获取每个符号（A,C,G,T,...）的总核苷酸计数、通过描述符字符串获取序列
  • 2021年3月2日：添加了“查找多个 FASTA 文件中的共同序列”的示例算法和新功能
© . All rights reserved.