如何拆分pacbio的数据

点击上方关注我们吧

随着pacbio测序的通量不断增加，样品的pooling也越来越常用。如何对混合测序的数据进行拆分呢，这是分析的第一步，也是很重要的一步。
lima是PacBio官方提供的拆分工具，在SMRT的v5.1.0版本中开始提供。替代之前版本的pbbarcode 和 bam2bam，可以为用户提供更好的使用体验。

支持的模式

支持每个样品有自己的barcode，pooling到一起，在一个cell上进行测序。支持以下四种不同的文库构建模式，获得的最终接头形式：

• 特异性的引物：通过一轮PCR，将barcode加在插入片段两端，然后连接上接头；

• barcode universal primer：通过两轮PCR，第一轮PCR将通用引物加在两端，之后再加barcode，最后做连接；

• barcode adapter：直接将barcode到序列两端；

• Probe-based linear barcoded adapters：通过连接酶直接链接上barcode和通用引物，然后用probe进行捕获，链接接头。

支持的barcode的顺序也各异，有以下三种形式都支持：

一条序列最多有两个barcode的区域，leading和tailing。可以只有单端barcode，也可以没有barcode。

对于对称的tailed文库设计，两端的barcode是一样的，不同的区别是barcode的方向。对于非对称的设计，两端的barcode是不一样的，因此需要提供两端的序列。

01

输入

一

输入文件

• 下机数据

  可以是CLR的序列，也可以用CCS的序列，bam格式。如果是RSII的数据，需要使用bax2bam进行转换成bam后，才能使用。

• barcode序列

  FASTA格式的文件。里面不能有重复，全部是大写字母，方向不敏感（正向和反向互补效果一样的），例如：

  >bc1000
  CTCTACTTACTTACTG
  >bc1001
  GTCGTATCATCATGTA
  >bc1002
  AATATACCTATCATTA

二

输入参数

• 文库构建相关的参数

• -s，--same ：两端barcode是一样的或者反向互补

• -d，--different：两端barcode不一样，只会输出--min-passes>=1的序列

• 数据类型相关参数

• --ccs：针对CCS后的数据使用，等价于-A 1 -B 4 -D 3 -I 3 -X 4 ，相当于增加了错配的罚分，并不是对数据进行ccs。

• --isoseq : Activate specialized IsoSeq mode。

• 二者不能同时使用

• 过滤参数

• --min-length ：拆分后，低于该长度的序列会被丢掉

• --max-input-lenght

• --min-score ：最低的barcode平均分数，默认是0，建议设置成26

• --min-end-score：最低的 end barcode的分数 ，在非对称barcode下有用。

• --min-ref-span 与 --min-scoring-region：前者是reference在barcode上区域的比例，默认是0.5；后者是最少的barcode region的个数；

• --min-passes

• 比对参数

• -A，-B，-D，-I，-X

• --peek：N，查看前N个ZMW，返回平均的barcode score。用于检测那些barcode正在使用。

• --guess：N，运行两次拆分。对此，所有的barcode都用于每个ZMW，计算每个barcode pair的出现次数，同时检测他们的barcode score是否大于N，只有这些barcode pair用于第二次拆分。在第二次拆分中，每个在第一轮中选择的pair用于拆分。与--peek结合起来使用。

• --guess-min-code：白名单中一个barcode需要的最少ZMW的个数。

• --peek-guess：默认是--peek 50000 --guess 45 --guess-min-count 10 ，这样可以提高运行效率。

02

输出

lima会输出多个文件，以相同的prefix作为输出文件的前缀。文件有：

• BAM：prefix.bam 包含去除后的序列，通过过滤标准的结果。

• Report：输出每个ZMW的比对信息，通过该文件可以追溯拆分过程，可惜header只能靠猜。一行是一个ZMW，如果使用--per-read，每一行是一个subread。

• summary：输出多少个ZMW被过滤掉，多少个接头是一致或者不一致，多少个序列。

• count：统计barcode对的个数和reads情况。

• Clipped region：输出被去掉的序列，使用--dump-clips。

• Removed records: 输出没有pass的序列或者没有接头的序列，使用--dump-removed。

• DataSet : 每个bam会有一个对应的说明文件。

• PBI : bam文件的索引文件。

03

质控工具

report_detail.R

对report文件进行可视化，并且进行质控

report_summary.R

对report文件，进行产量的绘图

04

运行模式

• 默认的情况下，lima使用ZMW对输入的数据进行分组，如果使用了--per-read，则会以每条序列为分组，进行分析。最终每个ZMW的结果是包含三个方面：对该ZMW中所有barcode区域的汇总，对提供的barcode序列中barcode对的结果；同一个ZMW的subread有相同的barcode和相同质量值。对于特定的barcode区域，每个barcode序列都会按照正向和反向互补进行比对，给出评分最高的结果。

• 如果两端是相同的barcode，那么使用--same 去掉不同的barcode对。

• 如果两端是不同的barcode，那么使用--different，去掉相同的barcode对，保留不同的。

05

评估表现

01

阳性预测率（PPV）

即拆分结果中，TP的比例是多少，反映了cross-calling的比例。为了计算PPV，不同长度和来源的扩增子，连接barcode后，测序 ，拆分，并且比对。通过这种方法，可以计算出TP和FP的个数。 许多影响拆分算法的因素也可以影响PPV的结果，例如barcode合成质量不高，不同孔的污染，建库的污染等。

同时，不同的barcode模式，双端是相同还是不同的barcode，barcode的种类的使用，都会ingxiangPPV。验证结果表明：

• barcode种类越多，PPV会下降

• 对称的barcode的PPV会好于非对称的文库结构

• 混合Sym和Asym会导致PPV更低，因此不建议采用

02

产量

lima去掉非期望的barcode对的结果，来提高PPV，这样会降低数据量。

非对称的产量会显著的低于对称的文库产量，因为对于非对称的文库，两端的序列都需要被测到，而对于对称的文库，只需要一端就可以了。

06

实现形式

使用SW算法，确定最佳的barcode得分和剪切位置，使用的是动态规划的方法。对于barcode的参考序列，使用的全局比对，而对于测序的序列使用局部比对，来确定得分的同时，确定clipping position。使用如下公式来确定barcode得分：

(100 * sw_score) / (sw_match_score * barcode_length)

完全匹配，是100分，得分越高，说明结果越准确。

结语

由于三代的单碱基的准确度不高，所以在拆分上不同于二代的算法，在效率上也慢了很多。期望未来应该有更快更好的软件出现。

1

reference

• https://github.com/PacificBiosciences/barcoding

作者：童蒙

编辑：amethyst