ubiome数据分析流程学习笔记1

uBiome 仅以 1%残值出售知识产权 2019 年 12 月 17 日，已经在 9 月 4 日向美国破产法院申请破产保护的 uBiome 公司宣布，将其名下 246 项微生物相关专利和 30 万人微生物组数据，以 770 万美元的价格转让给 Psomagen 公司，后者是韩国生物技术公司 Macrogen 在美国开展基因检测业务的子公司。

虽然ubiome倒下了，但是它的知识产权还是有价值的，要不也不会有公司出价收购，那么，我们看下它的方法具体细节是怎样的呢？

一般Next-Seq等illumina机器的最具性价比的测序读长是2＊150（150PE），而16S-V4区的长度是806R－515F＝291bp。那么，使用150PE测序方法对16S-V4区测序的话，获得的测序数据重叠区的长度是150＋150－291（包括引物长度）＝9bp，这个重叠理论上是可以拼接的，可是实际上，由于illumina测序反应（每次SBS会延伸一个碱基，大约耗时70分钟，据说之前拍照就需要40分钟。2013，https://blog.csdn.net/seallama/article/details/17613601） 推测大概在几十分钟，所以，150＊2个反应下来，试剂的衰减相当严重，测序最后获得的数据质量也急剧下降了，想要通过这最多9个的重叠碱基来实现拼接几乎是不可能完成的任务。

可是，这基本上是许多测序公司的测序方法，如果他们不采用特殊的方法是不能完成拼接，或者实现双向数据的有效利用的。其实这个问题，有两个方向可以解决，都已经发表了论文，一个是使用16S－V4区的通用引物来取代illumina的p5，p7测序引物，这样有效测序长度就变为300bp左右（不算上正反20bp左右的引物），重叠区域就变为50bp，许多论文和公司就是这么做的，这样做也比较简单和可靠，对于分析难度也比较小。如果一个lane全部测16S-V4，这是没问题的。

然而，16S-V4的测序数据量一般比较小，即使对于通量比较小的NextSeq－500机器，比如每次产出20GB算，一个样本获得测序rawdata 50Mb（量应该足够多了）算的话，至少需要500个样本混样测序（可能为了增加序列多样性，可能还要增加一定比例的phi噬菌体序列），而这对一般的检测机构，这么大的量是不可思议的。如果和别的不是测16S-V4文库一起测序的话就是不可行的，特别是16S－V4文库占比比较小的情况下。

这里，就引出了另一个方法，方法不好实现，就只有技术来扛了。另一个解决方案便是通过高水平的分析来解决不能拼接的问题了。我找到了几个方法，其中之一是读到的ubiome公司的一个方法。鉴于ubiome公司已经破产，据说是由于医保骗保，还有小道消息说是FBI盯上了遗传资源，但是从上面的报道来看，数据都可以卖给韩国人，小道消息应该是不正确。这个方法仅供参考，而且，由于他们公司注册了相关专利，也只能是学习下了，不可能商业应用。这篇文章发表在plos one，并不是多么高级的期刊，而且，这个公司的检测结果也面临质疑，所以，还是仅供参考。

数据分析流程

（其申请了专利，流程较复杂，特别是数据库要结合实验处理以减少假阳性和假阴性。）

1)数据库准备

a.首先从SILVA-16S数据库中找出能用V4通用引物扩出的序列，允许两个错配。检查得到的序列里是否有简并碱基，>20个以上的序列质量较差，由于采用双端测序，模拟测序产生的序列，去掉引物，就相当于正向有125bp和反向124bp接在一起。

b. 通过有选择地删除每个分类群的非特定扩增子集，可以为每个分类群创建几个经过筛选的数据库，并使用下面概述的步骤确定最佳数据库。采用100%的序列相似度和长度进行分析，排除不特异的扩增序列。然后，将注释到感兴趣分类群(Ti)和不同分类群(Dt)的序列按照dt/ti排列，并根据它们的商值分组。例如，在一个极端上，当商数为0时，我们创建了非常具体的数据库，其中不允许出现假阳性。在另一个极端，当商数很大时，例如，100，我们创建了非常敏感的数据库，在那里我们最大限度地识别真阳性，代价是产生假阳性。我们广泛地探索了这两个极端之间的不同可能的数据库，使用每个生成的数据库来预测每个感兴趣的物种和属的所有16S序列的分类。评估它们的性能，换句话说，选择性地从数据库中删除带有非特异性扩增的序列，同时最大限度地提高识别每个分类组中大多数序列的敏感性、特异性、精密度和阴性预测值。选择了灵敏度、特异度、精密度和阴性预测值均在90%以上、精密度和特异度之间的距离为最小可能值的分类单元作为每个分类单元的最佳数据库，目的是尽可能使精密度优于特异度。最终选择了28个种属进行分析。

2）物种注释

a.Q30过滤

b.去除引物和接头等

c. Swarm 2.15聚类，参数使用1个核苷酸的距离，fastidious和usearch-abundance。每个簇中最丰富的序列被认为是真正的生物序列，并被指定为簇中所有读取的计数。聚类中的其余读取被认为包含错误的测序的结果。

d. 使用vSearch从所有聚类中删除嵌合体。

e. 使用100%长度和100%同一性与来自手工处理过的silva数据库123版目标16Srna基因序列和分类注释的数据库进行比对。

f. 每个分类群的相对丰度是通过将与该分类群相关联的计数除以过滤后的总reads来确定的。

3.实验验证

从xTAG GPP (Luminex‘s xTAG Gastrointestinal Pathogen Panel)获得的样品, 对生物信息学流程进行了优化，以确保阳性结果真正意味着目标存在于样本中，并且仅在样本中没有目标存在时才获得阴性结果。

4. 健康队列的参考范围

为了确定28个种属的健康参考范围，我们建立了897健康个样本的队列。对于897个样本中的每一个，测定了微生物种群中每个目标的相对丰度。这些数据用来定义一个健康参考范围，例如艰难梭状芽孢杆菌在健康队列的~2%中发现，因此我们为其定义了一个相对丰度为0%到0.18%的健康范围。

5.不足之处

该方法不能识别肠球菌内不同的血清型，也不能检测能够区分致病性艰难梭菌或大肠杆菌与非致病性菌株的毒素基因，也不能在某些属水平的靶标中分辨物种。

6. 报告里的菌的情况

以上内容为文章中的，示例报告里的内容比文章中的有更新，菌有64属（种），35个种，29属。

阅读原文，分享给你这篇文章。https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0176555