GPU助力基因组重测序分析
作者：龚佳震
沐曦pde
背景介绍
1.1基因测序在癌症领域上的应用
癌症是目前人类所面临的最大敌人，其发病率随着年龄增长而升高，在人口老龄化加剧的当下，全社会的癌症负担也将愈发严峻。癌症难以治愈的原因之一是肿瘤具有强异质性，即相同症状、相同病理改变却可能由完全不同的基因变化而造成，以至于同类型癌症患者对于相同药物的药效反应有很大的差别 。这些基因上的变化包含了单碱基突变（snv）、小片段插入缺失（indel）、结构与拷贝数变异（sv & cnv）和甲基化等。随着近几年基因测序成本如图 1所示不断下降，在万元内即可完成人类的全基因组测序，gpu的技术发展也带来分析成本与时间的下降，于是用于检测基因组变化的重测序技术在癌症治疗中起到了越来越重要的作用。基因组重测序的应用主要有三个方向：
01基因测序技术提供肿瘤内基因分子水平的变化，为研究提供预见性，支持肿瘤新靶向药物研发。
02提高肿瘤早期诊断准确率，降低高危风险，从而提高癌症存活率。以在儿童中高发的儿童神经母细胞瘤（nb）为例，目前国内低、中危组nb患儿，长期无瘤生存率可达90%，但高危组nb患儿的5年总生存率不足30%。而如果能够做到早期诊断的话就可以通过手术切除联合化疗的方式进行治疗，从而提高治愈率。
03基于易感基因分子水平的基因变化检测，可以根据变化情况为患者提供个体化和预见性的治疗。
图 1 每兆数据dna测序开销
图中数据来源于wetterstrand ka. dna sequencing costs: data from the nhgri genome sequencing program (gsp) 和 next-generation-sequencing.v1.10.25
1.2基因测序手段介绍
当下基因组重测序常见的测序手段分为以illumina双端150bp为代表的短读长测序与以pacbio hifi或oxford nanopore technology（ont）为代表的长读长测序两种。如图1所示，illumina双端测序将dna打断两端加上不同接头与流动池中分布的接头结合不断合成、退火与扩增，扩增后利用带荧光的dntp与流动池中的dna链结合，由于不同的dntp发出不同荧光，于是可以通过扫描流动池中的荧光信号，得到序列信息；pacbio 将单分子与dna聚合酶固定在零模波导孔上，激光从孔的底部照入，激光不能穿过小孔且只能照亮孔底的一小部分区域。当dna聚合酶使用游离的dntp进行聚合反应时，dntp上的荧光基团被激光激发，通过传感器将光信号转化为电信号。ont则是通过让 dna 单分子通过带电的纳米孔，不同的碱基使得电流发生改变，传感器记录电流变化，测序仪将电流变化转化为碱基值。
图 2 测序原理示意图 (a) illumina 双端测序，
图片来源于dmlapato, illumina dye sequencing, https://en.wikipedia.org/w/index.php?title=illumina_dye_sequencing&oldid=1159305769#cite_note-illumina,_inc_2013-7
(b) pacbio 零模波导孔单分子实时测序，
图片来源于百迈客生物, pacbio三代测序仪, http://www.biomarker.com.cn/technology-services/pacbio
(c) ont 纳米孔测序，
图片来源于 kraft, long-read sequencing in human genetics, doi: 10.1007/s11825-019-0249-z
短读长测序的优势在于价格低廉与测序碱基准确度高，长读长测序拥有以下几个优点：
01长读长测序使用单分子测序，不需要扩增反应。一方面节省了时间，另一方面也避免了扩增中引入的错误，比如扩增效率不同带来的偏好性问题。
02长读长测序相比于短读长测序更容易跨过基因组的重复区域。在短读长测序结果中，至少748个与人类相关的基因存在测序暗黑区域，这个区域序列的read会出现低深度、低覆盖度及低比对质量。而造成这个现象的主要原因是在参考基因组中该区域的序列多次重复。而测序序列长到能够跨过重复区域的长读长测序正好解决了这个问题。
03长读长测序相比于短序列可以更容易得到核酸的甲基化水平。短读长测序基于荧光信号，通常需要将样本分为两份，一份利用亚硫酸氢盐处理将未甲基化的c转化为u，另一份直接测序进行比较获得基因组甲基化位点。而基于电信号的三代测序可以直接区分出甲基化的c与未甲基化的c，避免了亚硫酸氢盐等药品处理和多次测序引入的错误。
而近几年来测序的技术有两个明显的发展趋势：
01短读长测序和长读长测序成本均逐步降低。短读长测序成本降低幅度变缓，长读长测序的价格随着仪器的升级已经降到了三年前短读长测序成本。在一些应用场景中，比如结构变异的检测，由于长读长测序只需要三分之一短读长测序的数据量，长读长测序的成本已经低于短读长测序。
02随着试剂与算法的研发，长读长测序的准确度也逐渐接近短读长测序。在单碱基质量和测序价格的保证下，长序列能保留更原始的基因组特征。
于是目前研究所选择的测序手段逐渐从短读长测序向长读长测序转移。
1.3基因组重测序的流程介绍
如图3所示，重测序的流程主要包含以下大步骤：
图 3 基因组重测序流程与软件示意图。
图中紫色矩形部分为主要分析步骤，黄色矩形部分为分析类别，绿色矩形部分为软件，蓝色上标的软件为利用到了机器学习的软件。
01实验设计、分子实验与上机测序：这一步包括了需求分析和选择测序手段。分子实验通常包括核酸提取、核酸分子打断和样品的制备。上机测序将样品加入测序仪中，测序仪产生序列信息，完成测序。
02数据预处理：通常测序仪会利用传感器捕获电信号，测序仪的basecalling步骤和甲基化检测步骤将电信号转化为文本形式的序列。基因组测序的结果保存检测得到每一个核酸分子的碱基和该碱基的质量值，ont测序和pacbio测序也会保存原始电信号。甲基化检测的结果会同时输出每一个碱基的甲基化概率。长序列检测会发生随机的测序错误，被认为准确度不高，需要进行碱基修正，修正测序错误来提升碱基质量值。
03比对：这一步骤将测序得到的序列与参考基因组比对，确定该序列来自于参考基因组的哪个部分。rna测序得到的序列相比于dna缺失了内含子部分，长读长测序与短读长测序在错误率、复杂度上有区别。所以软件上会根据样本的类型和测序序列长度做区分。
04变异检测：变异检测将获得样本与参考基因组的差异，这些差异包括了单碱基的变异（snv）、小的插入缺失（indel）以及大的结构变异（sv）。短读长测序由于单碱基的准确度高，在snv和indel上有优势，而对于长度大于短读长测序序列长度的结构变异而言，长读长测序更有优势。所以软件上会根据需要检测的变异长度和测序序列长度做区分。
05后续分析：后续分析则包括了变异的过滤、注释、比较以及关联分析等过程，这一步将以万计的变异中筛选出与研究课题息息相关的变异，会跟随着课题不同而不同。
而这五个大步骤之中，第一步的时间改进则依赖于实验设计的改良以及试剂的研发。而后续步骤则需要依赖于算法的改进和计算设备的提升，在本文中将会介绍这些分析的常用软件和算法，并指出这些软件均可以通过对gpu和并行算法的优化进行提速。比如carpi等人利用合理实验设计和gpu，使得恶性疟原虫在计算成本减少4.6倍的同时整体分析速度提升27倍。
二测序数据预处理
2.1利用gpu获得碱基序列
将测序过程中产生的电信号转化为碱基的步骤被称为 basecalling。以 ont测序为例，如图 4所示，当核酸分子通过芯片上纳米孔时，生成不断变化的电信号，测序仪记录下电信号，并转化为碱基。
图 4 ont basecalling 过程示意图
图片来源于what is oxford nanopore technology (ont)sequencing?,https://www.yourgenome.org/facts/what-is-oxford-nanopore-technology-ont-sequencing/
由于碱基通过纳米孔的速度并不匀速，电流信号与预测得到的碱基单调但非对齐，为了解决这个问题，ont判断碱基的序列则是通过语音识别中常用的 ctc 作为损失函数来实现。ont开源basecalling 软件dorado 的r10.4.1系列模型中，模型分为编码器和解码器两个部分，解码器的模型如图 3所示，其接续线性条件随机场（crf）；解码器则是用viterbi算法进行的crf解码。dorado对每一个测序仪器构建fast、hac 和 sup 三种模型，模型依次增大且得到的碱基结果逐渐准确。r10.4.1模型中，在fast 模式下，解码器进入到线性crf层的特征向量为96，状态长度为3；hac模式为特征向量128，状态长度为4，而sup模式特征向量达到了256，状态长度为5。在 fast模型下gpu加速效果可以达到60倍以上。
图 5 dorado r10.4.1 fast模式的模型示意图
2.2deepconsensus深度学习修正测序错误
由于pacbio会有完全随机的错误发生（主要为短片段的插入与缺失），pacbio使用同一个零模波导孔内smrt 圈多次测序得到多条subreads序列的方法去除测序错误，矫正得到一个低错误率高质量值的序列，平均smrt 圈的测序圈数决定了测序的时间。pacbio revio测序仪使用ccs 和deepconsensus 模型来矫正测序错误获得高质量的序列。ccs基于隐马尔可夫模型，其gpu版本可加速至12倍。deepconsensus是一个编码器-only的transformer模型，与使用ccs的结果相比，deepconsensus能够得到更准确的多聚核苷酸位点。
deepconsensus将测序得到的reads切割，每份 100bp，每一份向量包含了序列、链方向、脉冲宽度、脉冲间隔以及信噪比等信息。通过一个6个自注意力层的编码器后，使用以softmax作为激活函数的前馈神经网络解码得到真实的序列。由于gpu在神经网络训练和推理上的优势，相比于cpu版本的deepconsensus，gpu加速效果能够达到 3.3倍以上。
图 6 deepconsensus 流程示意图
图片来源于baid et al. deepconsensus improves the accuracy of sequences with a gap-aware sequence transformer. doi: 10.1038/s41587-022-01435-7
三测序序列比对
两条序列比对问题是序列比对的基础，以图7所示，序列1 gactttcc 与 序列2 tagactacc 进行比对得到比对表。两条dna序列的比对通常有四种状态：插入、缺失、匹配和错配，而匹配和错配往往会在比对问题中合并为同一个状态。
图 7 seq1与seq2的序列比对示意图
smith-waterman算法和 needleman-wunsch算法是序列比对中常用的算法，其中smith-waterman算法用于局部比对而needleman-wunsch算法应用于全局序列比对。而实际案例中，由于参考基因组的序列长，比如人的参考序列总长度在3gb，一次测序深度10x测序所得reads数目在2百万左右，直接使用局部比对的smith-waterman算法在时间成本上是不可接受的，于是有两种加速方案：第一种为各测序所得序列间并行运算，另一种则是对现有算法进行改进，minimap2等基needleman-wunsch的比对软件便应运而生。在本节中首先介绍smith-waterman与needleman-wunsch的动态规划算法与其gpu加速情况，然后介绍minimap2特殊的gpu加速。
3.1smith-waterman与needleman-wunsch 算法与gpu加速
两个算法拥有相同的步骤：首先初始化得分矩阵，根据相似性矩阵和罚分矩阵填满表格，后根据分数回溯得到比对结果。如图 8所示，gotoh提出的比对方式具体如下：
01确定罚分矩阵，罚分矩阵定义了碱基匹配的奖励分s、开启间隔区域（插入或者缺失）的惩罚分q 和持续间隔e 的惩罚分。这一步往往根据测序的特性进行确定。
02根据罚分矩阵来初始化打分矩阵，通常会将分数存储在三个打分矩阵 h、e 和f 中，分别代表着匹配分数、插入分数和缺失分数。smith-waterman算法中，打分矩阵h、e 和f 的首列与首行初始化为0；而needleman-wunsch算法中，打分矩阵f首列初始化公式为f0,j=f0,j-1-e，打分矩阵e首行初始化公式为ei,0=ei-1,0-e，打分矩阵h的首列和首行初始化为hi,j=max{ei,j,fi,j}
03计算打分矩阵h算法中的每一个位置的分数hi,j 。smith-waterman 算法中hi,j=max{hi-1,j-1+s,ei,j,fi,j,0} ，而needleman-wunsch算法中，hi,j=max{hi-1,j-1+s,ei,j,fi,j}，其中 ei,j=max{hi-1,j-q,ei-1,j}-e,fi,j=max{hi,j-1-q,fi,j-1}-e。
当打分矩阵h、e 和f 全部填写完成后，需要回溯来获取最优的比对序列。smith-waterman算法从打分矩阵中的最大值开始回溯直到hi,j=0 为止；而needleman-wunsch算法从打分矩阵h 右下角开始回溯，直到左上角为止。回溯路径即为比对路径。此时对于长度m 和长度n 的序列对而言，时间复杂度为o(mn)。
myers提出了bit-vector加速该算法，打分矩阵hi,j 中的依赖于hi-1,j-1、hi-1,j 和hi,j-1，如果以列向量的角度看的话，如图 9所示，那么hj 仅依赖于hj-1。该算法将计算打分矩阵hi,j 变为计算行差值δhi,j 与列差值δvi,j。将时间复杂度降为o(n)。
图 9 bit-vector 加速示意图
siriwardena等人提出了另一个加速方案，如图 10所示，次对角线上的各个元素没有相互依赖，故拓展到次对角线的block没有依赖，因此通过各个block之间的并行计算进行gpu加速。该方案运行时间与比对序列长度成非线性正相关，在如人染色体的数量级上能够达到40倍的加速效果。
图 10 次对角线分块加速示意图
与cpu 和 gpu执行时间比较图
除了上述提到的方案以外，使用动态规划算法的时候需要频繁用到max、min等计算，gpu同时采用特殊的硬件指令来加速动态规划算法，预计可以加速7倍。
3.2比对软件的新星minimap2的流程与gpu加速
为了解决smith-waterman局部比对在大规模数据中耗时长的问题，li开发了minimap2 软件极大减少了比对的时间。minimap2 是在比对领域较新且广泛使用的软件，相比于其它软件如bwa-mem、star等有特定使用场景的限制，minimap2的使用范围比较广泛，其适用于二代的dna短读长测序数据以及三代dna长读长测序数据，也可以被应用rna-seq这种有大间隔的数据。如图11所示，minimap2将测序得到的 reads与参考基因组的比对的过程主要分为 seeding、chaining 和 alignment 三个步骤。
图 11 minimap2 流程示意图
01seeding这一步快速定位参考基因组与read之间完全匹配，无插入与缺失的短序列，得到若干有匹配位置信息的锚点（anchors）。锚点使用a(x,y,w)表示，意味着参考基因组的坐标[x-w+1, x]与read的坐标[y-w+1, y]完全匹配，每一条read会得到若干个锚点。
02chaining这一步则是将seeding 这步得到的锚点连起来，确定这条read在参考基因组上的真实位置。在 minimap2中，chaining是用一个一维的动态规划算法计算f 解决的，将所有锚点根据参考基因组的位置坐标x进行排序，第 i个anchor的分数f(i)可以由公式 f(i)=max{max{f(j)+α(j,i)-β(j,i),wi}(i>j≥1)} , 计算得到，其中 α(j,i)=min{min{yi-yj,xi-xj},wi} 作为奖励分。为两个锚点之间相同的长度；β(j,i)=γc ((yi-yj )-(xi-xj ))，作为两个锚点之间发生插入或者缺失间隔的罚分，γc 是一个与间隔长度相关的函数，可以调整 γc 以适应不同测序类型。minimap2也会设置一个最长间隔g，当max{yi-yj,xi-xj}>g时，罚分β(j,i)=∞。若有n个锚点，那么直接计算最大f时间复杂度为o(n2)，基于大部分f(i) 的最大值在i附近，故只会计算h 次，于是minimap2将时间控制在o(hn)内。对于i每次得到f 后，会回溯到 j 直到f(j)=wj ，此时将回溯过程中的锚点记为“已用”，得到一条链，每一个链根据匹配长度和间隔长度计算比对分数s。当我们将所有的锚点标记为“已用”后，可以得到所有链和其分数s，分数最高的链便是得到的最优链。
图 12 minimap2不同步骤所花时间比例图。图片来源于sadasivan et al. accelerating minimap2 for accurate long read alignment on gpus. doi: 10.26502/jbb.2642-
03alignment 这一步，则是将测序得到的 read 与 chaining 这一步得到的链序列进行 base to base 的全局比对。这一步则可以用上一小节所提到的needleman-wunsch算法完成。
minimap2的耗时时间如图 12所示，chaining和 alignment 是耗时最长的时间。因此提升 chaining 和 alignment步骤的速度是提升 miniamp2 运行速度的关键。
sadasivan 等人提出了minimap2 chaining步骤的并行方案。如图 13 所示，他们首先将 minimap2 的chaining步骤的逆向搜索变成顺向搜索，并且对于每一个i 锚点并行计算其顺向的 >i 锚点是否为最大分数。gpu通过这种方式并行即可以将minimap2提速3.8倍左右。
图 13 mm2-ax 加速 minimap2 chaining 的方案
图片来源于sadasivan et al. accelerating minimap2 for accurate long read alignment on gpus. doi: 10.26502/jbb.2642-91280067
四变异检测
4.1snv与indel的检测
在生物学上，通常需要关注的是测序得到的变异。单碱基变异（snv）是常见的一种变异方式，snv 是单碱基的变异，碱基的突变有时候会改变氨基酸，或者使得转录提前终止。snv的检测主要分为两个步骤：1. 确定并过滤变异位置； 2. 确定基因型。在本节中将介绍两种软件，一种是以 haplotypecaller和mutect2为代表的基于传统概率学模型建模的软件，另一种则是以deepvariant为代表的利用深度学习进行检测的软件。
4.1.1概率模型 haplotypecaller / mutect2检测流程
haplotypecaller 与 mutect2 是 gatk 中重要的 snv 与 indel calling 的软件。两款软件在一定程度上拥有相同的原理，在概率模型上的区别，导致haplotypecaller 主要用于性细胞变异的检测，而mutect2 更适用于体细胞变异的检测。如图 14两款软件都在确定变异区间后，利用局部组装得到的单倍型并将测序得到的reads用smith-waterman算法比对到单倍型中来减少测序、比对等错误信息，通过 pairhmm 来得到每一个单倍型的似然值，然后根据各自的概率模型获得最终的基因型。
图 14 haplotypecaller与mutect2 的通用流程示意图
每一个单倍型的似然值时haplotypecaller和mutect2的关键步骤，即求单倍型的集合h 对于所有测序得到的reads的集合r的比对结果集合a的概率，即 p(r│h)=∑a p(r,a|h)，那么只需要求解每一个比对的概率即可。其中我们已经看到了重新比对中smith-waterman的gpu加速方案，haplotypecaller中的pairhmm这一步亦可以用gpu并行加速。
pairhmm 的过程与前述的比对算法相似，定义了三个状态: 匹配m, 插入i 和缺失d。如图15所示，三个状态会发生转移, m→i 或者 m→d，即从匹配状态转化为插入状态或者缺失状态概率为δ、 i→m 从插入状态结束间隔变回匹配状态概率为ε、d→m从缺失状态结束间隔变回匹配状态概率为ε，也会发生持续 m→m 即延续匹配状态概率为1-2δ、i→i 延续插入状态概率为1-ε 和d→d 延续缺失状态，概率为1-ε。基于生物的序列特征，往往1-ε50bp的序列，与snv或小片段的indel相比，一旦变异发生通常更不可逆。结构变异的检测有着重要的意义：
1由于不可逆的特性，结构变异往往是人种区别、物种分化的关键。
2稀有的结构变异往往与疾病或者癌症发生关联。
3在肿瘤细胞中变异会出现大量的高密度的异常的结构变异，甚至会出现染色体碎裂的现象。
结构变异的检测的方法如图 19所示，包括了基于双端测序方向的read pair、基于测序深度变化的 read depth、基于比对时异常比对的split reads 和将测序得到的 reads组装与参考基因组比对的 assembly。目前短序列有manta、delly，而长序列有 sniffles2、cutesv等软件，基于比对时异常比对的split reads是主流的方法。
图 19 结构变异检测的常用方法
图片来源于alkan et al. genome structural variation discovery and genotyping. doi: 10.1038/nrg2958
4.2.1svision深度学习模型
如 sniffles2、manta等结构变异软件，过滤得到比对异常的reads，将比对异常点根据其在参考基因组的位置、异常的长度等方式进行聚类，每一个类即获得一个结构变异位点。但是这缺少了发现复合结构变异的能力。通常复合结构变异检测都依赖于现有的模式，如易位后缺失等，缺乏一个统一的结构。
svision 便是为了解决这个问题而出现的模型31。如图 20所示，svision包含了编码器、tmor和解释器。编码器将测序的序列转化为图片，其识别支持变异序列和参考基因组序列的碱基通过var-ref和ref-ref图像，从var-ref中去除ref-ref中已有的重复序列，以提高后续分析的准确度。在tmor步骤中，svison将包含多个复合结构变异断点的图像通过5层卷积层、1层展平层和2层全连接层的卷积神经网络进行预测发生变异的概率。最后，解释器描述了不同的复合结构变异。由于这是一个卷积神经网络的模型，gpu的矩阵乘法上的优势能够比cpu版本的svision快1.8倍。
图 20 svision 流程示意图
图片来源于 lin et al. svision: a deep learning approach to resolve complex structural variants. doi: 10.1038/s41592-022-01609-w
四总结与展望
变异检测是定位功能基因或疾病基因的基础。如表格 1所示，从测序数据的产生、预处理、基因组的比对、变异检测都可以使用gpu加速来提升检测的速度。高速分析能够节省测序的时间成本，从而降低测序价格和分析成本。这可以让研究人员更愿意进行测序来建立与完善疾病肿瘤相关的数据库，为将来攻克疾病与肿瘤打下基础。
表格 1 文中所提的软件gpu与cpu运行时间比
可以看到，在目前基因测序领域，传统软件算法因为设备局限性，在面对复杂的生物系统时，选择使用简化的概率模型建模来牺牲准确度而获得更快的速度。而随着gpu算力的进步，越来越多研究人员选择深度学习或者大模型来得到更精确的结果。除了文中所提到的数据预处理、基因组比对和变异检测步骤以外，在后续分析中，机器学习也起到了更重要的作用，
比如用线性模型建模关联样本表型与基因型的gwas可以使用 deepnull 得到更好的关联结果。
我们预计，随着机器学习算法和gpu设备的发展，快速而准确的基因测序可以帮助研发人员设计靶向药物，帮助普通民众更早发现潜在风险，帮助医生更准确地对病情进行诊断治疗，精准医疗未来可期。
名词列表
贸泽电子与运动控制公司Trinamic 签署全球分销协议
betBEB钱包正在让加密数字货币走向现实
PLC与计算机的数据通信
智能公交与车载视频一体化终端的研究与展望
华为mate10/mate10pro即将发布:倪飞亲自下场宣传,确认有三个版本,保时捷版买不起,全面屏版售价良心
GPU助力基因组重测序分析
华为确定无法和谷歌合作了即将启动备选方案
基地继电器的原理图
多款互联网存款产品下架，是否对银行有影响？
iPhone X导入OLED面板，韩系、日系、陆系厂商在OLED面板布局
奥迪将增加未来五年电动汽车预算
2021年第十六届“中国芯”集成电路产业促进大会 暨“中国芯”优秀产品征集结果发布仪式在珠海开幕
NVIDIA公布一条“大新闻：NVIDIA DLSS支持深度学习抗锯齿技术
我国液压、气压动力机械企业营收呈增长态势，未来发展趋势如何
与Arm抢市场，国内这家嵌入式CPU IP核供应商冲刺科创板
瑞丰恒S9型紫外激光器的应用优势及助力仪表盘打标
英伟达Win10 GeForce 461.09驱动可能会导致崩溃
人形机器人索菲亚，正以量产之势进入人类社会
一加5什么时候上市?一加手机5真机图现身:金属机身+骁龙835+首款8GB手机
过载过流的基础知识汇总
沐曦pde
背景介绍
1.1基因测序在癌症领域上的应用
癌症是目前人类所面临的最大敌人，其发病率随着年龄增长而升高，在人口老龄化加剧的当下，全社会的癌症负担也将愈发严峻。癌症难以治愈的原因之一是肿瘤具有强异质性，即相同症状、相同病理改变却可能由完全不同的基因变化而造成，以至于同类型癌症患者对于相同药物的药效反应有很大的差别 。这些基因上的变化包含了单碱基突变（snv）、小片段插入缺失（indel）、结构与拷贝数变异（sv & cnv）和甲基化等。随着近几年基因测序成本如图 1所示不断下降，在万元内即可完成人类的全基因组测序，gpu的技术发展也带来分析成本与时间的下降，于是用于检测基因组变化的重测序技术在癌症治疗中起到了越来越重要的作用。基因组重测序的应用主要有三个方向：
01基因测序技术提供肿瘤内基因分子水平的变化，为研究提供预见性，支持肿瘤新靶向药物研发。
02提高肿瘤早期诊断准确率，降低高危风险，从而提高癌症存活率。以在儿童中高发的儿童神经母细胞瘤（nb）为例，目前国内低、中危组nb患儿，长期无瘤生存率可达90%，但高危组nb患儿的5年总生存率不足30%。而如果能够做到早期诊断的话就可以通过手术切除联合化疗的方式进行治疗，从而提高治愈率。
03基于易感基因分子水平的基因变化检测，可以根据变化情况为患者提供个体化和预见性的治疗。
图 1 每兆数据dna测序开销
图中数据来源于wetterstrand ka. dna sequencing costs: data from the nhgri genome sequencing program (gsp) 和 next-generation-sequencing.v1.10.25
1.2基因测序手段介绍
当下基因组重测序常见的测序手段分为以illumina双端150bp为代表的短读长测序与以pacbio hifi或oxford nanopore technology（ont）为代表的长读长测序两种。如图1所示，illumina双端测序将dna打断两端加上不同接头与流动池中分布的接头结合不断合成、退火与扩增，扩增后利用带荧光的dntp与流动池中的dna链结合，由于不同的dntp发出不同荧光，于是可以通过扫描流动池中的荧光信号，得到序列信息；pacbio 将单分子与dna聚合酶固定在零模波导孔上，激光从孔的底部照入，激光不能穿过小孔且只能照亮孔底的一小部分区域。当dna聚合酶使用游离的dntp进行聚合反应时，dntp上的荧光基团被激光激发，通过传感器将光信号转化为电信号。ont则是通过让 dna 单分子通过带电的纳米孔，不同的碱基使得电流发生改变，传感器记录电流变化，测序仪将电流变化转化为碱基值。
图 2 测序原理示意图 (a) illumina 双端测序，
图片来源于dmlapato, illumina dye sequencing, https://en.wikipedia.org/w/index.php?title=illumina_dye_sequencing&oldid=1159305769#cite_note-illumina,_inc_2013-7
(b) pacbio 零模波导孔单分子实时测序，
图片来源于百迈客生物, pacbio三代测序仪, http://www.biomarker.com.cn/technology-services/pacbio
(c) ont 纳米孔测序，
图片来源于 kraft, long-read sequencing in human genetics, doi: 10.1007/s11825-019-0249-z
短读长测序的优势在于价格低廉与测序碱基准确度高，长读长测序拥有以下几个优点：
01长读长测序使用单分子测序，不需要扩增反应。一方面节省了时间，另一方面也避免了扩增中引入的错误，比如扩增效率不同带来的偏好性问题。
02长读长测序相比于短读长测序更容易跨过基因组的重复区域。在短读长测序结果中，至少748个与人类相关的基因存在测序暗黑区域，这个区域序列的read会出现低深度、低覆盖度及低比对质量。而造成这个现象的主要原因是在参考基因组中该区域的序列多次重复。而测序序列长到能够跨过重复区域的长读长测序正好解决了这个问题。
03长读长测序相比于短序列可以更容易得到核酸的甲基化水平。短读长测序基于荧光信号，通常需要将样本分为两份，一份利用亚硫酸氢盐处理将未甲基化的c转化为u，另一份直接测序进行比较获得基因组甲基化位点。而基于电信号的三代测序可以直接区分出甲基化的c与未甲基化的c，避免了亚硫酸氢盐等药品处理和多次测序引入的错误。
而近几年来测序的技术有两个明显的发展趋势：
01短读长测序和长读长测序成本均逐步降低。短读长测序成本降低幅度变缓，长读长测序的价格随着仪器的升级已经降到了三年前短读长测序成本。在一些应用场景中，比如结构变异的检测，由于长读长测序只需要三分之一短读长测序的数据量，长读长测序的成本已经低于短读长测序。
02随着试剂与算法的研发，长读长测序的准确度也逐渐接近短读长测序。在单碱基质量和测序价格的保证下，长序列能保留更原始的基因组特征。
于是目前研究所选择的测序手段逐渐从短读长测序向长读长测序转移。
1.3基因组重测序的流程介绍
如图3所示，重测序的流程主要包含以下大步骤：
图 3 基因组重测序流程与软件示意图。
图中紫色矩形部分为主要分析步骤，黄色矩形部分为分析类别，绿色矩形部分为软件，蓝色上标的软件为利用到了机器学习的软件。
01实验设计、分子实验与上机测序：这一步包括了需求分析和选择测序手段。分子实验通常包括核酸提取、核酸分子打断和样品的制备。上机测序将样品加入测序仪中，测序仪产生序列信息，完成测序。
02数据预处理：通常测序仪会利用传感器捕获电信号，测序仪的basecalling步骤和甲基化检测步骤将电信号转化为文本形式的序列。基因组测序的结果保存检测得到每一个核酸分子的碱基和该碱基的质量值，ont测序和pacbio测序也会保存原始电信号。甲基化检测的结果会同时输出每一个碱基的甲基化概率。长序列检测会发生随机的测序错误，被认为准确度不高，需要进行碱基修正，修正测序错误来提升碱基质量值。
03比对：这一步骤将测序得到的序列与参考基因组比对，确定该序列来自于参考基因组的哪个部分。rna测序得到的序列相比于dna缺失了内含子部分，长读长测序与短读长测序在错误率、复杂度上有区别。所以软件上会根据样本的类型和测序序列长度做区分。
04变异检测：变异检测将获得样本与参考基因组的差异，这些差异包括了单碱基的变异（snv）、小的插入缺失（indel）以及大的结构变异（sv）。短读长测序由于单碱基的准确度高，在snv和indel上有优势，而对于长度大于短读长测序序列长度的结构变异而言，长读长测序更有优势。所以软件上会根据需要检测的变异长度和测序序列长度做区分。
05后续分析：后续分析则包括了变异的过滤、注释、比较以及关联分析等过程，这一步将以万计的变异中筛选出与研究课题息息相关的变异，会跟随着课题不同而不同。
而这五个大步骤之中，第一步的时间改进则依赖于实验设计的改良以及试剂的研发。而后续步骤则需要依赖于算法的改进和计算设备的提升，在本文中将会介绍这些分析的常用软件和算法，并指出这些软件均可以通过对gpu和并行算法的优化进行提速。比如carpi等人利用合理实验设计和gpu，使得恶性疟原虫在计算成本减少4.6倍的同时整体分析速度提升27倍。
二测序数据预处理
2.1利用gpu获得碱基序列
将测序过程中产生的电信号转化为碱基的步骤被称为 basecalling。以 ont测序为例，如图 4所示，当核酸分子通过芯片上纳米孔时，生成不断变化的电信号，测序仪记录下电信号，并转化为碱基。
图 4 ont basecalling 过程示意图
图片来源于what is oxford nanopore technology (ont)sequencing?,https://www.yourgenome.org/facts/what-is-oxford-nanopore-technology-ont-sequencing/
由于碱基通过纳米孔的速度并不匀速，电流信号与预测得到的碱基单调但非对齐，为了解决这个问题，ont判断碱基的序列则是通过语音识别中常用的 ctc 作为损失函数来实现。ont开源basecalling 软件dorado 的r10.4.1系列模型中，模型分为编码器和解码器两个部分，解码器的模型如图 3所示，其接续线性条件随机场（crf）；解码器则是用viterbi算法进行的crf解码。dorado对每一个测序仪器构建fast、hac 和 sup 三种模型，模型依次增大且得到的碱基结果逐渐准确。r10.4.1模型中，在fast 模式下，解码器进入到线性crf层的特征向量为96，状态长度为3；hac模式为特征向量128，状态长度为4，而sup模式特征向量达到了256，状态长度为5。在 fast模型下gpu加速效果可以达到60倍以上。
图 5 dorado r10.4.1 fast模式的模型示意图
2.2deepconsensus深度学习修正测序错误
由于pacbio会有完全随机的错误发生（主要为短片段的插入与缺失），pacbio使用同一个零模波导孔内smrt 圈多次测序得到多条subreads序列的方法去除测序错误，矫正得到一个低错误率高质量值的序列，平均smrt 圈的测序圈数决定了测序的时间。pacbio revio测序仪使用ccs 和deepconsensus 模型来矫正测序错误获得高质量的序列。ccs基于隐马尔可夫模型，其gpu版本可加速至12倍。deepconsensus是一个编码器-only的transformer模型，与使用ccs的结果相比，deepconsensus能够得到更准确的多聚核苷酸位点。
deepconsensus将测序得到的reads切割，每份 100bp，每一份向量包含了序列、链方向、脉冲宽度、脉冲间隔以及信噪比等信息。通过一个6个自注意力层的编码器后，使用以softmax作为激活函数的前馈神经网络解码得到真实的序列。由于gpu在神经网络训练和推理上的优势，相比于cpu版本的deepconsensus，gpu加速效果能够达到 3.3倍以上。
图 6 deepconsensus 流程示意图
图片来源于baid et al. deepconsensus improves the accuracy of sequences with a gap-aware sequence transformer. doi: 10.1038/s41587-022-01435-7
三测序序列比对
两条序列比对问题是序列比对的基础，以图7所示，序列1 gactttcc 与 序列2 tagactacc 进行比对得到比对表。两条dna序列的比对通常有四种状态：插入、缺失、匹配和错配，而匹配和错配往往会在比对问题中合并为同一个状态。
图 7 seq1与seq2的序列比对示意图
smith-waterman算法和 needleman-wunsch算法是序列比对中常用的算法，其中smith-waterman算法用于局部比对而needleman-wunsch算法应用于全局序列比对。而实际案例中，由于参考基因组的序列长，比如人的参考序列总长度在3gb，一次测序深度10x测序所得reads数目在2百万左右，直接使用局部比对的smith-waterman算法在时间成本上是不可接受的，于是有两种加速方案：第一种为各测序所得序列间并行运算，另一种则是对现有算法进行改进，minimap2等基needleman-wunsch的比对软件便应运而生。在本节中首先介绍smith-waterman与needleman-wunsch的动态规划算法与其gpu加速情况，然后介绍minimap2特殊的gpu加速。
3.1smith-waterman与needleman-wunsch 算法与gpu加速
两个算法拥有相同的步骤：首先初始化得分矩阵，根据相似性矩阵和罚分矩阵填满表格，后根据分数回溯得到比对结果。如图 8所示，gotoh提出的比对方式具体如下：
01确定罚分矩阵，罚分矩阵定义了碱基匹配的奖励分s、开启间隔区域（插入或者缺失）的惩罚分q 和持续间隔e 的惩罚分。这一步往往根据测序的特性进行确定。
02根据罚分矩阵来初始化打分矩阵，通常会将分数存储在三个打分矩阵 h、e 和f 中，分别代表着匹配分数、插入分数和缺失分数。smith-waterman算法中，打分矩阵h、e 和f 的首列与首行初始化为0；而needleman-wunsch算法中，打分矩阵f首列初始化公式为f0,j=f0,j-1-e，打分矩阵e首行初始化公式为ei,0=ei-1,0-e，打分矩阵h的首列和首行初始化为hi,j=max{ei,j,fi,j}
03计算打分矩阵h算法中的每一个位置的分数hi,j 。smith-waterman 算法中hi,j=max{hi-1,j-1+s,ei,j,fi,j,0} ，而needleman-wunsch算法中，hi,j=max{hi-1,j-1+s,ei,j,fi,j}，其中 ei,j=max{hi-1,j-q,ei-1,j}-e,fi,j=max{hi,j-1-q,fi,j-1}-e。
当打分矩阵h、e 和f 全部填写完成后，需要回溯来获取最优的比对序列。smith-waterman算法从打分矩阵中的最大值开始回溯直到hi,j=0 为止；而needleman-wunsch算法从打分矩阵h 右下角开始回溯，直到左上角为止。回溯路径即为比对路径。此时对于长度m 和长度n 的序列对而言，时间复杂度为o(mn)。
myers提出了bit-vector加速该算法，打分矩阵hi,j 中的依赖于hi-1,j-1、hi-1,j 和hi,j-1，如果以列向量的角度看的话，如图 9所示，那么hj 仅依赖于hj-1。该算法将计算打分矩阵hi,j 变为计算行差值δhi,j 与列差值δvi,j。将时间复杂度降为o(n)。
图 9 bit-vector 加速示意图
siriwardena等人提出了另一个加速方案，如图 10所示，次对角线上的各个元素没有相互依赖，故拓展到次对角线的block没有依赖，因此通过各个block之间的并行计算进行gpu加速。该方案运行时间与比对序列长度成非线性正相关，在如人染色体的数量级上能够达到40倍的加速效果。
图 10 次对角线分块加速示意图
与cpu 和 gpu执行时间比较图
除了上述提到的方案以外，使用动态规划算法的时候需要频繁用到max、min等计算，gpu同时采用特殊的硬件指令来加速动态规划算法，预计可以加速7倍。
3.2比对软件的新星minimap2的流程与gpu加速
为了解决smith-waterman局部比对在大规模数据中耗时长的问题，li开发了minimap2 软件极大减少了比对的时间。minimap2 是在比对领域较新且广泛使用的软件，相比于其它软件如bwa-mem、star等有特定使用场景的限制，minimap2的使用范围比较广泛，其适用于二代的dna短读长测序数据以及三代dna长读长测序数据，也可以被应用rna-seq这种有大间隔的数据。如图11所示，minimap2将测序得到的 reads与参考基因组的比对的过程主要分为 seeding、chaining 和 alignment 三个步骤。
图 11 minimap2 流程示意图
01seeding这一步快速定位参考基因组与read之间完全匹配，无插入与缺失的短序列，得到若干有匹配位置信息的锚点（anchors）。锚点使用a(x,y,w)表示，意味着参考基因组的坐标[x-w+1, x]与read的坐标[y-w+1, y]完全匹配，每一条read会得到若干个锚点。
02chaining这一步则是将seeding 这步得到的锚点连起来，确定这条read在参考基因组上的真实位置。在 minimap2中，chaining是用一个一维的动态规划算法计算f 解决的，将所有锚点根据参考基因组的位置坐标x进行排序，第 i个anchor的分数f(i)可以由公式 f(i)=max{max{f(j)+α(j,i)-β(j,i),wi}(i>j≥1)} , 计算得到，其中 α(j,i)=min{min{yi-yj,xi-xj},wi} 作为奖励分。为两个锚点之间相同的长度；β(j,i)=γc ((yi-yj )-(xi-xj ))，作为两个锚点之间发生插入或者缺失间隔的罚分，γc 是一个与间隔长度相关的函数，可以调整 γc 以适应不同测序类型。minimap2也会设置一个最长间隔g，当max{yi-yj,xi-xj}>g时，罚分β(j,i)=∞。若有n个锚点，那么直接计算最大f时间复杂度为o(n2)，基于大部分f(i) 的最大值在i附近，故只会计算h 次，于是minimap2将时间控制在o(hn)内。对于i每次得到f 后，会回溯到 j 直到f(j)=wj ，此时将回溯过程中的锚点记为“已用”，得到一条链，每一个链根据匹配长度和间隔长度计算比对分数s。当我们将所有的锚点标记为“已用”后，可以得到所有链和其分数s，分数最高的链便是得到的最优链。
图 12 minimap2不同步骤所花时间比例图。图片来源于sadasivan et al. accelerating minimap2 for accurate long read alignment on gpus. doi: 10.26502/jbb.2642-
03alignment 这一步，则是将测序得到的 read 与 chaining 这一步得到的链序列进行 base to base 的全局比对。这一步则可以用上一小节所提到的needleman-wunsch算法完成。
minimap2的耗时时间如图 12所示，chaining和 alignment 是耗时最长的时间。因此提升 chaining 和 alignment步骤的速度是提升 miniamp2 运行速度的关键。
sadasivan 等人提出了minimap2 chaining步骤的并行方案。如图 13 所示，他们首先将 minimap2 的chaining步骤的逆向搜索变成顺向搜索，并且对于每一个i 锚点并行计算其顺向的 >i 锚点是否为最大分数。gpu通过这种方式并行即可以将minimap2提速3.8倍左右。
图 13 mm2-ax 加速 minimap2 chaining 的方案
图片来源于sadasivan et al. accelerating minimap2 for accurate long read alignment on gpus. doi: 10.26502/jbb.2642-91280067
四变异检测
4.1snv与indel的检测
在生物学上，通常需要关注的是测序得到的变异。单碱基变异（snv）是常见的一种变异方式，snv 是单碱基的变异，碱基的突变有时候会改变氨基酸，或者使得转录提前终止。snv的检测主要分为两个步骤：1. 确定并过滤变异位置； 2. 确定基因型。在本节中将介绍两种软件，一种是以 haplotypecaller和mutect2为代表的基于传统概率学模型建模的软件，另一种则是以deepvariant为代表的利用深度学习进行检测的软件。
4.1.1概率模型 haplotypecaller / mutect2检测流程
haplotypecaller 与 mutect2 是 gatk 中重要的 snv 与 indel calling 的软件。两款软件在一定程度上拥有相同的原理，在概率模型上的区别，导致haplotypecaller 主要用于性细胞变异的检测，而mutect2 更适用于体细胞变异的检测。如图 14两款软件都在确定变异区间后，利用局部组装得到的单倍型并将测序得到的reads用smith-waterman算法比对到单倍型中来减少测序、比对等错误信息，通过 pairhmm 来得到每一个单倍型的似然值，然后根据各自的概率模型获得最终的基因型。
图 14 haplotypecaller与mutect2 的通用流程示意图
每一个单倍型的似然值时haplotypecaller和mutect2的关键步骤，即求单倍型的集合h 对于所有测序得到的reads的集合r的比对结果集合a的概率，即 p(r│h)=∑a p(r,a|h)，那么只需要求解每一个比对的概率即可。其中我们已经看到了重新比对中smith-waterman的gpu加速方案，haplotypecaller中的pairhmm这一步亦可以用gpu并行加速。
pairhmm 的过程与前述的比对算法相似，定义了三个状态: 匹配m, 插入i 和缺失d。如图15所示，三个状态会发生转移, m→i 或者 m→d，即从匹配状态转化为插入状态或者缺失状态概率为δ、 i→m 从插入状态结束间隔变回匹配状态概率为ε、d→m从缺失状态结束间隔变回匹配状态概率为ε，也会发生持续 m→m 即延续匹配状态概率为1-2δ、i→i 延续插入状态概率为1-ε 和d→d 延续缺失状态，概率为1-ε。基于生物的序列特征，往往1-ε50bp的序列，与snv或小片段的indel相比，一旦变异发生通常更不可逆。结构变异的检测有着重要的意义：
1由于不可逆的特性，结构变异往往是人种区别、物种分化的关键。
2稀有的结构变异往往与疾病或者癌症发生关联。
3在肿瘤细胞中变异会出现大量的高密度的异常的结构变异，甚至会出现染色体碎裂的现象。
结构变异的检测的方法如图 19所示，包括了基于双端测序方向的read pair、基于测序深度变化的 read depth、基于比对时异常比对的split reads 和将测序得到的 reads组装与参考基因组比对的 assembly。目前短序列有manta、delly，而长序列有 sniffles2、cutesv等软件，基于比对时异常比对的split reads是主流的方法。
图 19 结构变异检测的常用方法
图片来源于alkan et al. genome structural variation discovery and genotyping. doi: 10.1038/nrg2958
4.2.1svision深度学习模型
如 sniffles2、manta等结构变异软件，过滤得到比对异常的reads，将比对异常点根据其在参考基因组的位置、异常的长度等方式进行聚类，每一个类即获得一个结构变异位点。但是这缺少了发现复合结构变异的能力。通常复合结构变异检测都依赖于现有的模式，如易位后缺失等，缺乏一个统一的结构。
svision 便是为了解决这个问题而出现的模型31。如图 20所示，svision包含了编码器、tmor和解释器。编码器将测序的序列转化为图片，其识别支持变异序列和参考基因组序列的碱基通过var-ref和ref-ref图像，从var-ref中去除ref-ref中已有的重复序列，以提高后续分析的准确度。在tmor步骤中，svison将包含多个复合结构变异断点的图像通过5层卷积层、1层展平层和2层全连接层的卷积神经网络进行预测发生变异的概率。最后，解释器描述了不同的复合结构变异。由于这是一个卷积神经网络的模型，gpu的矩阵乘法上的优势能够比cpu版本的svision快1.8倍。
图 20 svision 流程示意图
图片来源于 lin et al. svision: a deep learning approach to resolve complex structural variants. doi: 10.1038/s41592-022-01609-w
四总结与展望
变异检测是定位功能基因或疾病基因的基础。如表格 1所示，从测序数据的产生、预处理、基因组的比对、变异检测都可以使用gpu加速来提升检测的速度。高速分析能够节省测序的时间成本，从而降低测序价格和分析成本。这可以让研究人员更愿意进行测序来建立与完善疾病肿瘤相关的数据库，为将来攻克疾病与肿瘤打下基础。
表格 1 文中所提的软件gpu与cpu运行时间比
可以看到，在目前基因测序领域，传统软件算法因为设备局限性，在面对复杂的生物系统时，选择使用简化的概率模型建模来牺牲准确度而获得更快的速度。而随着gpu算力的进步，越来越多研究人员选择深度学习或者大模型来得到更精确的结果。除了文中所提到的数据预处理、基因组比对和变异检测步骤以外，在后续分析中，机器学习也起到了更重要的作用，
比如用线性模型建模关联样本表型与基因型的gwas可以使用 deepnull 得到更好的关联结果。
我们预计，随着机器学习算法和gpu设备的发展，快速而准确的基因测序可以帮助研发人员设计靶向药物，帮助普通民众更早发现潜在风险，帮助医生更准确地对病情进行诊断治疗，精准医疗未来可期。
名词列表
贸泽电子与运动控制公司Trinamic 签署全球分销协议
betBEB钱包正在让加密数字货币走向现实
PLC与计算机的数据通信
智能公交与车载视频一体化终端的研究与展望
华为mate10/mate10pro即将发布:倪飞亲自下场宣传,确认有三个版本,保时捷版买不起,全面屏版售价良心
GPU助力基因组重测序分析
华为确定无法和谷歌合作了即将启动备选方案
基地继电器的原理图
多款互联网存款产品下架，是否对银行有影响？
iPhone X导入OLED面板，韩系、日系、陆系厂商在OLED面板布局
奥迪将增加未来五年电动汽车预算
2021年第十六届“中国芯”集成电路产业促进大会 暨“中国芯”优秀产品征集结果发布仪式在珠海开幕
NVIDIA公布一条“大新闻：NVIDIA DLSS支持深度学习抗锯齿技术
我国液压、气压动力机械企业营收呈增长态势，未来发展趋势如何
与Arm抢市场，国内这家嵌入式CPU IP核供应商冲刺科创板
瑞丰恒S9型紫外激光器的应用优势及助力仪表盘打标
英伟达Win10 GeForce 461.09驱动可能会导致崩溃
人形机器人索菲亚，正以量产之势进入人类社会
一加5什么时候上市?一加手机5真机图现身:金属机身+骁龙835+首款8GB手机
过载过流的基础知识汇总