1. 课题背景#

5-甲基胞嘧啶（m5C）作为一种重要的表观遗传标记，不仅存在于DNA中，在RNA中也扮演着关键角色。这种化学修饰对基因表达调控机制产生深远影响，影响细胞功能、发育过程和疾病状态。随着高通量测序技术的快速发展，我们现在拥有了前所未有的分辨率和规模来探索m5C的分布模式及其生物学意义。然而，尽管高通量测序为m5C检测提供了有力支持，其应用并非没有挑战。

传统的m5C检测方法，如RNA-Bis-seq和UBS-seq，虽然能够以单碱基分辨率识别m5C位点，但通常依赖于C→T转换策略来实现。这种方法的一个主要局限性是容易导致高假阳性率，因为一些正常的C到T变异可能被错误解释为m5C的存在。因此，在降低假阳性的同时提高m5C检测的准确性和可靠性已成为该领域的关键挑战之一。

本任务旨在优化一个系统的工作流程，该流程整合了一套专门用于检测高通量测序数据中RNA m5C的生物信息学软件包，目标是实现高效、精确的m5C位点识别。具体来说，该工作流程从原始测序数据（FASTQ文件）开始，依次执行接头修剪、低质量序列过滤、rRNA/tRNA过滤、基因组比对、m5C位点检测和严格的数据过滤，最终生成高置信度的m5C候选位点列表。通过这种方法，我们不仅可以有效降低假阳性率，还能显著提高计算效率，从而为后续的功能验证和机制探索提供可靠的数据基础。通过系统地优化整个m5C检测流程，我们可以更深入地理解m5C在转录后调控中的作用，揭示其在健康和疾病状态下的复杂动态变化，并为未来的研究奠定坚实的基础。这对于推动精准医学、疾病诊断和治疗具有重要意义。在以下章节中，我们将详细解释对此工作流程所做的优化。

2. 集群配置#

我们在三台计算设备上进行了优化工作，每台设备的硬件配置相同，如下表所示。

以下是我们的软件环境配置：

首先，使用 GCC-12.2.0 作为基准编译器，并按默认配置编译软件。通过执行推荐的脚本并利用 time 命令精确测量各阶段耗时，我们获得了数据处理流程（阶段1）的性能数据：case1 耗时 8,476 秒，case2 耗时 11,146 秒，case3 耗时 9,944 秒。鉴于实验平台使用 Intel® Xeon® Gold 6348 处理器，我们切换到 Intel oneAPI 系列编译器，以充分利用硬件优势。这一选择不仅能利用 Intel 编译器对其架构的高度优化特性，还能为高级矢量扩展指令集提供有力支持，显著提升计算密集型任务的执行效率。此外，oneAPI 工具包中的 Vtune Profiler 为代码性能分析和瓶颈识别提供了强大工具。

将编译器更换为 Intel(R) oneAPI DPC++/C++ Compiler 2022.2.1 后，我们对软件配置进行了相应优化：在 CFLAGS 和 CPPFLAGS 中添加了优化选项 “-O3 -mavx2 -funroll-loops -qopenmp”；为 samtools 的依赖库 htslib 添加了 “—with-libdeflate” 编译选项以增强压缩性能；并将测试脚本中的 gzip 替换为效率更高的 libdeflate-gzip。经过这些优化，阶段1的处理时间显著减少：case1 降至 8,002 秒，case2 降至 10,636 秒，case3 降至 9,421 秒。

3. 流程分析与初步尝试#

为方便起见，我们将 SRR23538290、SRR23538291 和 SRR23538292 对应的文件分别称为 case1、case2 和 case3，并在下文中使用这些名称。

我们对阶段1中主要工具的耗时进行了详细提取和量化分析，并根据获得的数据绘制了相应的饼图，如图 1 所示。从图 1 可以清楚地观察到，整个数据处理流程中耗时最长的三个关键步骤是：原始测序数据与双参考序列的比对、基于唯一分子标识符的去重，以及 3n 转换表的生成。其中，序列比对过程主要依赖 hisat-3n 工具，它通过高效的索引结构和比对算法实现大规模测序数据的快速比对。UMI 去重操作使用 UMICollapse（一个 Java 包）实现，它能有效识别和去除 PCR 扩增过程中引入的重复序列，从而提高数据准确性。3n 转换表的生成是进一步处理比对结果的关键步骤。基于以上分析，hisat-3n 和 UMICollapse 的性能直接影响整个工作流程的效率。因此，优化这两个工具的性能将成为未来研究工作的重点。

图1：阶段1各步骤耗时占比

考虑到我们使用的 CPU 核心数（56）远超推荐脚本中设置的核心数，我们最初将脚本中支持并行计算的步骤的核心参数调整为 56。然而，经过相当长的执行时间后，我们才最终获得结果。随后，我们严格按照推荐脚本的默认配置重新运行实验，发现使用 56 核心的脚本执行时间明显高于推荐脚本。这一现象促使我们深刻反思软件的并行化效率，也为我们后续的性能分析和优化工作提供了重要背景和动力。

4. 文件缓存#

在推荐脚本中，阶段1运行时间最长的操作是一个组合命令，涉及 Samtools、hisat-3n-table、cut 和 gzip 等工具，如下面的代码所示。这些工具通过管道在彼此之间传输数据，一个工具的输出直接作为下一个工具的输入。具体来说，该组合命令通过顺序连接工具实现高效的数据流和处理。然而，这种串行数据处理方式可能成为性能瓶颈，尤其是在处理大规模测序数据时。管道数据传输的开销和工具间的依赖关系可能导致整体效率下降。因此，优化此组合命令，特别是探索其并行化潜力和提高资源利用率，将是未来性能优化的关键重点之一。

1
(time samtools view -e "rlen<100000" -h $OUTPUT/7/SRR23538290.mRNA.genome.mapped.sorted.dedup.bam | hisat-3n-table -p 16 -u --alignments - --ref /path/to/Homo_sapiens.GRCh38.dna.primary_assembly.fa --output-name /dev/stdout --base-change C,T | cut -f 1,2,3,5,7 | libdeflate-gzip -c > $OUTPUT/9/SRR23538290_unfiltered_uniq.tsv.gz) > $LOG/9/log.txt 2>&1

为了单独评估 hisat-3n-table 的性能，我们移除了原有的管道机制，改为将每个步骤的输出写入磁盘文件。后续步骤再读取这些生成的文件进行数据处理。修改后的代码如下：

1
{ time samtools view -e "rlen<100000" -h $OUTPUT/7/SRR23538290.mRNA.genome.mapped.sorted.dedup.bam > $OUTPUT/9/SRR23538290.filtered.bam; } 2> $LOG/9/log.txt
2
{ time hisat-3n-table -p 16 -u --alignments $OUTPUT/9/SRR23538290.filtered.bam --ref /path/to/Homo_sapiens.GRCh38.dna.primary_assembly.fa --output-name /dev/stdout --base-change C,T > $OUTPUT/9/SRR23538290.hisat_output.tsv; } 2>> $LOG/9/log.txt
3
{ time cut -f 1,2,3,5,7 $OUTPUT/9/SRR23538290.hisat_output.tsv > $OUTPUT/9/SRR23538290.filtered_columns.tsv; } 2>> $LOG/9/log.txt
4
{ time libdeflate-gzip -c $OUTPUT/9/SRR23538290.filtered_columns.tsv > $OUTPUT/9/SRR23538290_unfiltered_uniq.tsv.gz; } 2>> $LOG/9/log.txt

实验结果出乎意料：命令拆分后的运行时间比原始脚本缩短了近一半。以 case1 为例，该命令的执行时间从 1,751 秒显著下降到 990 秒。我们推测性能提升的主要原因是原始管道机制传输的数据量过大，导致 I/O 瓶颈。将中间结果写入磁盘有效缓解了这个问题，因为文件读写操作能更好地利用系统的缓冲机制，从而减少数据传输开销。通过分析中间文件大小，我们进一步验证了这一假设的正确性——中间文件确实占用了大量磁盘空间，这也部分解释了为什么在管道机制下性能受到限制。这一发现为未来的优化工作提供了宝贵指导，强调了设计数据流过程以最小化 I/O 开销、提高整体工作效率的重要性。

5. hisat-3n-table 模块优化#

1
void append(int threadID) {
2
    string* line;
3
    Alignment newAlignment;
4

5
    while (working) {
6
        workerLock[threadID]->lock();
7
        if(!linePool.popFront(line)) {
8
            workerLock[threadID]->unlock();
9
            this_thread::sleep_for (std::chrono::nanoseconds(1));
10
            continue;
11
        }
12
        while (refPositions.empty()) {
13
            this_thread::sleep_for (std::chrono::microseconds(1));
14
        }
15
        newAlignment.parse(line);
16
        returnLine(line);
17
        appendPositions(newAlignment);
18
        workerLock[threadID]->unlock();
19
    }
20
}

1
while (alignmentFile->good()) {
2
    positions->getFreeStringPointer(line);
3
    if (!getline(*alignmentFile, *line)) {
4
        positions->returnLine(line);
5
        break;
6
    }
7

8
    if (line->empty() || line->front() == '@') {
9
        positions->returnLine(line);
10
        continue;
11
    }
12
    // limit the linePool size to save memory
13
    while(positions->linePool.size() > 1000 * nThreads) {
14
        this_thread::sleep_for (std::chrono::microseconds(1));
15
    }
16
    // if the SAM line is empty or unmapped, get the next SAM line.
17
    if (!getSAMChromosomePos(line, samChromosome, samPos)) {
18
        positions->returnLine(line);
19
        continue;
20
    }
21
    // if the samChromosome is different than current positions' chromosome, finish all SAM line.
22
    // then load a new reference chromosome.
23
    if (samChromosome != positions->chromosome) {
24
        // wait all line is processed
25
        while (!positions->linePool.empty() || positions->outputPositionPool.size() > 100000) {
26
            this_thread::sleep_for (std::chrono::microseconds(1));
27
        }
28
        positions->appendingFinished();
29
        positions->moveAllToOutput();
30
        positions->loadNewChromosome(samChromosome);
31
        reloadPos = loadingBlockSize;
32
        lastPos = 0;
33
    }
34
    // if the samPos is larger than reloadPos, load 1 loadingBlockSize bp in from reference.
35
    while (samPos > reloadPos) {
36
        while (!positions->linePool.empty() || positions->outputPositionPool.size() > 100000) {
37
            this_thread::sleep_for (std::chrono::microseconds(1));
38
        }
39
        positions->appendingFinished();
40
        positions->moveBlockToOutput();
41
        positions->loadMore();
42
        reloadPos += loadingBlockSize;
43
    }
44
    if (lastPos > samPos) {
45
        cerr << "The input alignment file is not sorted. Please use sorted SAM file as alignment file." << endl;
46
        throw 1;
47
    }
48
    positions->linePool.push(line);
49
    lastPos = samPos;
50
}

1
while (alignmentFile->good())
2
{
3
    if (!getline(*alignmentFile, line_1))
4
    {
5
        break;
6
    }
7

8
    if (line_1.empty() || line_1.front() == '@' || !getSAMChromosomePos(line_1, samChromosome, samPos))
9
    {
10
        continue;
11
    }
12

13
    V_SamPos.emplace_back(samPos);
14

15
    if (samChromosome != chromosome)
16
    {
17
        V_Index.emplace_back(count);
18
        V_SamChromosome.emplace_back(samChromosome);
19
        chromosome = samChromosome;
20
    }
21

22
    linePool.emplace_back(line_1);
23
    count++;
24
}

1
for(long long int i_index = 0; i_index < V_Index.size() - 1; i_index++)
2
{
3
    positions->loadNewChromosome(V_SamChromosome[i_index]);
4

5
    auto temp_reloadPos = loadingBlockSize;
6

7
    for(long long int j_index = V_Index[i_index]; j_index < V_Index[i_index + 1]; ++j_index)
8
    {
9
        while (V_SamPos[j_index] > temp_reloadPos)
10
        {
11
            positions->moveBlockToOutput();
12
            positions->loadMore();
13
            temp_reloadPos += loadingBlockSize;
14
        }
15
        positions->append();
16
    }
17
    positions->moveAllToOutput();
18
}

1
*out_ << "ref\tpos\tstrand\tconvertedBaseQualities\tconvertedBaseCount\tunconvertedBaseQualities\tunconvertedBaseCount\n";
2
Position* pos;
3
for(int i = 0; i < V_Index.size() - 1; ++i)
4
{
5
    while (outputPositionPool[i].popFront(pos))
6
    {
7
        *out_ << pos->chromosome << '\t'
8
                << to_string(pos->location) << '\t'
9
                << pos->strand << '\t'
10
                << pos->convertedQualities << '\t'
11
                << to_string(pos->convertedQualities.size()) << '\t'
12
                << pos->unconvertedQualities << '\t'
13
                << to_string(pos->unconvertedQualities.size()) << '\n';
14
    }
15
}

1
char* b_1;
2
for (int i = 0; i < line_2.size(); i++) {
3
    Position newPos_1;
4
    newPos_1.set(targetChromosome, location + i);
5
    b_1 = &line_2[i];
6
    if (CG_only) {
7
        if (lastBase == 'C' && *b_1 == 'G') {
8
            refPositions.back().set('+');
9
            newPos_1.set('-');
10
        }
11
    } else {
12
        if (*b_1 == convertFrom) {
13
            newPos_1.set('+');
14
        } else if (*b_1 == convertFromComplement) {
15
            newPos_1.set('-');
16
        }
17
    }
18
    refPositions.emplace_back(newPos_1);
19
    lastBase = *b_1;
20
}
21

22
/* 原始代码
23
Position* newPos_1;
24
char* b_1;
25
for (int i = 0; i < line_2.size(); i++) {
26
    if (!freePositionPool.popFront(newPos_1)) {
27
        newPos_1 = new Position();
28
    }
29
    newPos_1->set(targetChromosome, location + i);
30
    b_1 = &line_2[i];
31
    if (CG_only) {
32
        if (lastBase == 'C' && *b_1 == 'G') {
33
            refPositions.back()->set('+');
34
            newPos_1->set('-');
35
        }
36
    } else {
37
        if (*b_1 == convertFrom) {
38
            newPos_1->set('+');
39
        } else if (*b_1 == convertFromComplement) {
40
            newPos_1->set('-');
41
        }
42
    }
43
    refPositions.emplace_back(newPos_1);
44
    lastBase = *b_1;
45
}
46
*/

6. 数据复用#

在流程优化过程中，实验观察发现，连续的 hisat-3n-table 模块都将来自 samtools 的相同预处理数据作为其输入源。为了解决这个冗余操作，我们实现了一个中间文件缓存机制来优化资源分配。具体来说，通过将 samtools 的初始处理结果持久存储为中间文件，我们确保一次数据过滤操作就能满足所有后续相关工作流程的输入需求。该策略有效消除了重复执行资源密集型过滤操作的需要，将 I/O 操作从两次减少到一次。结果，系统资源消耗显著降低。实验验证证实，这种利用中间文件复用的优化方法不仅保持了数据处理的一致性，还将整体工作流程效率提高了约 23%，同时减少了约 37% 的临时文件存储空间。

7. 并发操作#

如前所述，我们对 hisat-3n-table 模块进行了优化，显著提高了其多核性能。然而，随着线程数逐渐增加到 16，进一步增加线程数带来的额外改进微乎其微，表明程序已经以极快的速度运行。为了充分利用 56 核处理器的性能，我们决定并发执行四个 hisat-3n-table 命令，并并行运行两个 samtools 操作，因为它们没有数据依赖关系。

具体来说，我们在每个命令的末尾附加一个 & 符号，使其能够在后台运行而不阻塞后续步骤。最后，我们使用 wait 命令阻塞，直到所有后台进程完成。脚本流程概述如下：

1
# echo "------------11-------------"
2
(time samtools view -@ 8 -e "[XM] * 20 <= (qlen-sclen) && [Zf] <= 3 && 3 * [Zf] <= [Zf] + [Yf]" ...) > $LOG/11/log.txt 2>&1 &
3

4
echo "------------9-------------"
5
{ time samtools view -e "rlen<100000" ...; } 2> $LOG/9/log.txt &
6

7
wait
8

9
{ time samtools view -e "rlen<100000" ...; } 2> $LOG/12/log.txt
10
{ time hisat-3n-table ...; } 2>> $LOG/9/log.txt &
11
{ time hisat-3n-table ...; } 2>> $LOG/10/log.txt &
12
{ time hisat-3n-table ...; } 2>> $LOG/12/log.txt &
13
{ time hisat-3n-table ...; } 2>> $LOG/13/log.txt &
14

15
wait

8. UMI 去重优化#

当使用 UMICollapse 工具对成功比对到人类基因组的 BAM 文件进行 UMI 去重时，我们观察到 CPU 资源利用率显著，只有一个核心持续达到高占用率。不幸的是，该工具利用多核并行计算的能力并未得到有效发挥。最初尝试使用 umitools 作为替代方案发现，它不仅缺乏并行计算选项，而且运行效率明显低于 UMICollapse。因此，去重操作继续依赖 UMICollapse 作为处理工具。

基于人类基因组参考序列的染色体分布特征，我们设计并实现了一种基于染色体坐标的并行优化策略。该策略的实施涉及以下关键步骤：

1. 染色体信息提取： 利用 samtools idxstats 命令，获取整个基因组中所有染色体的名称，为后续基于染色体的文件拆分提供基础数据支持。
2. 特殊染色体处理： 识别包括 X、Y 染色体和线粒体 DNA（非 1-22 号染色体）在内的特殊染色体，并将其名称存储在一个列表（special_chrom_list）中。鉴于这些特殊染色体的比对记录数量远少于标准染色体，我们使用 samtools view 命令从源文件中提取并合并这些记录，以提高处理效率。
3. 标准染色体数据拆分： 对于比对到 1-22 号染色体的测序数据，这些数据通常涉及更大的数据集，我们实施了一种基于染色体编号独立拆分数据的策略。每条染色体的比对记录都存储在单独的文件中。实验结果表明，使用 samtools view 进行拆分所需的时间可以忽略不计（< 总处理时间的 2%），显示出高效率。
4. 并行处理： 数据拆分后，我们采用并行计算策略，同时启动多个 UMICollapse 实例，对每个特定染色体的文件进行 UMI 去重。
5. 结果整合： 最后，我们使用 samtools merge 命令整合并重建所有处理后的结果。值得注意的是，由于排序已在早期阶段完成，此合并步骤不会打乱文件顺序。通过多次实验验证，这种拆分合并策略在确保结果准确性的同时，显著提高了处理效率。

此步骤的具体工作流程如图 4 所示。

图4：步骤执行工作流程图

该优化方案的完整实现细节在 umi-java.sh 自动化脚本中提供。实验结果表明，引入多线程并行计算框架显著提高了系统资源利用率，并有效利用了多核能力，实现了约 6 倍的计算效率提升（基于 n=5 次独立实验）。在项目时间有限的情况下，与在工具源代码层面进行并行化修改相比，这种基于任务的并行化策略展示了卓越的工程实施效率，以及在资源利用率和实用性之间更好的平衡。

9. cutseq 和 hisat2-align-s 的协作并行#

该优化方案的完整实现细节在 umi-java.sh 自动化脚本中提供。实验结果表明，引入多线程并行计算框架显著提高了系统资源利用率，并有效利用了多核能力，实现了约 6 倍的计算效率提升（基于 n=5 次独立实验）。在项目时间有限的情况下，与在工具源代码层面进行并行化修改相比，这种基于任务的并行化策略展示了卓越的工程实施效率，以及在资源利用率和实用性之间更好的平衡。

1
/path/to/hisat2-align-s --wrapper basic-0 --index /path/to/Homo_sapiens.GRCh38.ncrna.fa --summary-file /path/to/map2ncrna.output.summary --new-summary -q -p 16 --base-change C,T --mp 8,2 --no-spliced-alignment --directional-mapping -U /path/to/SRR23538290.fastq_cut --3N

使用火焰图进行的性能分析表明，计算程序的核心热点函数是 align()，作为基因序列比对的顶层接口，初步确定为系统性能的潜在瓶颈。进一步研究发现，hybridSearch_recur() 函数消耗时间最多，它通过多模式搜索策略实现碱基对匹配，但受限于其复杂的递归机制。理论上，将递归范式转换为迭代模型可以减少堆栈内存使用和复杂性，从而提高执行效率。然而，实验数据表明，这种迭代重构并未显著提升性能。因此，我们将重点转移到系统性地优化次要热点，如 getAnchorHits()。详细检查其算法逻辑和实现，发现了诸如位操作和数据压缩技术等优化。通过多维度的优化实验，证实当前实现已接近最优性能，无法获得进一步显著提升。

图5：火焰图

为了解决当前计算密集型模块的优化瓶颈，我们重新关注 hisat-3n-table 流程中的 I/O 性能问题。使用 GDB，我们逐行跟踪分析了文件读写和数据转换过程。我们发现 readOrigBuf 缓冲区处理来自异常巨大的 FASTQ 格式拆分文件（fastq_cut，>10GB）的原始文本数据，这显著超过了典型文件大小。进一步分析执行流程发现，在 go() 函数的最外层控制循环中，存在对 fastq_cut 每一行的顺序读取操作（每次循环迭代一行）。尽管此读取操作仅用于信息输出，不贡献于核心计算，但其数十亿次的迭代使文件读取成为关键性能瓶颈。

计时分析实验表明，上游 cutseq 模块的生成速率（输出带宽）远超下游 hisat2-align-s 模块的消耗速率（输入带宽），导致显著的吞吐量不平衡。基于这些发现，我们提出一个基于生产者-消费者模型的并发执行框架：通过实现两阶段流水线并行化策略，确保 fastq_cut 文件的同步生成和消耗（设置 5 秒等待间隔，远高于实际 1-2 秒的文件创建延迟），从而有效重叠计算资源。以下是此优化方法的简化实现框架。

1
time cutseq ... &
2
sleep 5
3
time hisat-3n ... | samtools ... &
4
wait

通过应用上述方法，我们能够实现资源的优化分配，有效减轻初始 cutseq 步骤带来的运行时开销。然而，第四步 hisat2-align-s 的执行时间并未见显著改善。为了解决这个问题，我们深入分析了与 —index 参数对应的文件在程序中的使用方式。我们的发现表明，这些文件在初始化阶段被预加载到内存中，并在调用时加载到相关变量中。读取和调用机制设计合理，几乎没有明显的优化空间。

10. 总结#

三个案例的加速趋势几乎相同；因此，我们以 case2 为例说明各阶段的加速比，如图 6 所示。

图6：Case2 各阶段耗时与加速比

对于各个工具，最耗时的组件 hisat-3n-table 得到了显著优化，UMICollapse 的使用也有了实质性改进，这两个成为整个流程中加速贡献最高的部分。关于整体工作流程，对管道的修改和并发执行的引入不可或缺。结果，完成所有案例的阶段1和阶段2现在仅需 1773 秒，如图 7 所示。

图7：阶段1和阶段2总耗时