被特殊物种序列虚晃一枪的日子

以下内容，纯属虚构，如有雷同，哇哈哈哈。

前言

二代测序实属物种鉴定的利器，呼吸介入诊断体系将其纳入之后如虎添翼。除了操作前的临床理化影像分析，操作中的快速现场细胞学评价（Rapid On-Site Evaluation，ROSE）和操作后的宏基因组二代测序（Metagenomic Next Generation Sequencing，mNGS）大大提高了呼吸介入操作的诊断效率。但这些技术在应用过程中仍需要存在很多具体问题，需要在实践中不断探索解决。尤其是在发现疑似特殊类型病原体时，需要严谨的验证，并在这个过程中不断完善分析流程。

缘起

病房接连入住两名患者，均因为基础病情而处于免疫抑制状态，肺内病变较为复杂，鉴别诊断当中包括了病毒性感染不能除外。一名患者病情进展非常迅猛，另一名患者病情已经进入吸收恢复期。按照我的操作常规，灌洗刷检活检一应俱全，各路标本分别送检，其中就包括病变部位的支气管肺泡灌洗液（Broncho-Alveolar Lavage Fluid，BALF）送检mNGS。检测公司非常给力，出具检测报告的同时还会提供下机数据，客户支持做得非常到位，不像其他有些公司店大欺客，张口闭口要这要那，完全不给自己提供任何业务机会。

由于患者病情加重，各级查房都非常重视。拿到数据之后，就按照之前《宏基因组单个样本数据处理流程笔记》当中所述流程进行分析。最后点开花冠图，就觉得出现了一个值得注意的物种，占据病毒界的比例还不低，加之临床和影像表现需要将病毒感染列入鉴别诊断范围，这就引起了注意。

回溯

当前面临一个问题：这是否就是由该病原体引起的感染？

在前面一篇文章《从kraken2注释结果中提取指定物种序列》当中，我就分享了如何将指定物种序列从CleanData当中提取出来的圡办法。将这些提取的序列再次验证，这次使用最传统的blast。使用的时候才发现blast还能针对特定物种进行匹配，我要指定的物种信息还被放在了首页，非常方便。除此之外，kraken2的开发者还在github上提供了一些特殊物种的序列供研究交流。于是我下载了这些序列，使用本地的blast进行比较。使用makeblastdb命令对上述序列进行建库，blastn命令进行序列比较，所有提取序列对于特殊物种的一致性高达95%以上。

为了怕自己的验证有误，于是先将这些提取的序列送到一个专门做数据处理的师妹那里，她验证结果跟我一样。后来我还将原始的数据发给她从头验证，不知道现在是否已经得出验证结果。

不管怎样，既然存在这样的序列，而且一致性还那么高，我也第一时间向上级汇报。

辨析

第三名患者的测序结果也跟前两人相同。情况已经上报，至于如何应对，肉食者谋，而对于前一个问题的探索仍然继续进行。

分布

并不是说只要序列能够匹配到基因组上就算该病原体真的存在，在基因组上的分布情况同样重要。既然二代测序在提取核酸之后是将其随机片段化，那么理论上来说，物种片段匹配到基因组上时应该分布比较均匀。于是我从blast上观察了一下这些片段在目标物种基因组上的分布情况。这些序列分布非常集中，分布于整个29+k基因组的28+k的N蛋白编码和21+k的S蛋白编码上，位置特别集中，只有偶尔一两条序列分布在前面13+k的ORF1ab。

为了观察其他病毒序列分布的特点，专门把一名巨细胞病毒肺炎患者的序列按照上面的方法给提取出来。对这些序列在CMV基因组上进行匹配后，发现这些位置相对比较分散，并不会特别集中于某一两点。但是考虑到该特殊病原体本身的特性相对CMV这种DNA病毒来说更容易降解，所以还特地找了一个身家相似的鼻病毒（Rhinovirus）来比较，直接观察到的结论也和CMV相同。

这些数据是采用76bp读长的Illumina体系进行检测，深度较为有限。为了担心测序深度不够导致没能充分覆盖基因组，我请公司使用另一个50bp读长的更高通量体系复查上述患者的RNA。结果……序列竟然消失了，而其他物种基本还存在。但这些序列的一致性却非常高，所以有必要检查一下来源。

核酸

对于微生物样本的检测分为DNA和RNA两种。鼻病毒属于RNA病毒，但在DNA检测的序列中同样存在。毕竟，鼻病毒除了自己直接合成蛋白之外，要进行核酸复制就需要逆转录到DNA再合成新的RNA链。而巨细胞病毒这样的DNA病毒，也会产生RNA来合成蛋白。所以无论是DNA检测还是RNA检测，都应该出现这些病原体的踪迹。但是，这些特殊物种的序列并没有出现在DNA检测中，却仅仅出现在RNA检测中，而且相对丰度非常高，这不能不让人怀疑其来源是否是真的病原体。

搜捕

为了确证我的猜想，在一个时机恰当的夜晚，我集中分析了截止目前所有的测序数据。我一直使用Intel的NUC第十代寒霜峡谷来做基因数据处理，速度特别慢，直接起名叫“亀”。所以数据量一大就犯懒，毕竟键盘一敲这12条线程就不知道跑到啥时候去了。没办法，因为太穷，脸都要不起了。这次难得豁出去了，这老多的数据，通宵也得干完，得把问题搞清楚。

到凌晨三点多，三十多名患者共几十G的数据终于处理完毕。唉，也就是干活儿前讲课讲得脑子有点儿抽抽，干到两点多钟才想起为啥自己要一遍一遍的倒腾，写个Shell不香么……算了，剩下一两个了已经。今年使用76bp读长体系检测的数据就存在上述问题，即DNA检测正常，而RNA检测中就出现了特殊物种。而且这些序列片段都集中于上述两个区域，并没有在基因组的其他部位出现。

至此可以推测，序列并非来源于标本，而是来源于RNA的反应体系，最大可能是来自试剂。

真相

确认上述情况之后，集中上述提取的核酸序列，跟公司技术进行沟通。技术表示已经经过技术手段确定这段序列来自于试剂的某个特定组分，并非完整的病原基因组，而只有其中极短的一部分，这也是为何总共只检测到两个位点总长度不足200bp的原因。而对于信息学处理部分，只需要将这一段序列进行比对后予以去除即可，无论是bwa还是bowtie还是别的工具，都可以。

并非终章

特殊物种的序列在面前虚晃一枪，还好虚惊一场。不过从中能获得的技术经验确实很宝贵，特别是对于物种注释之后再次鉴定确认，还需要进行逆向提取之后对于目标微生物的基因组进行分析，判断结果的真实性。技术流程还需要不断迭代完善，才能从容应对复杂的临床状况。这些流程和相关分析代码我会陆续整理之后与大家分享交流。

最后感谢大家听我讲了这么啰嗦的一个故事，总结陈词如下：

以上内容，纯属虚构，如有雷同，哇哈哈哈。