1. 技术产生背景

PacBio构建这种文库的一个主要原因是，在其测序芯片中有很多零模波导孔（ZMW），每个孔在测序时只能测一条DNA分子，这样其测序通量就受到孔数和文库长度的限制。近些年PacBio确实也在不停的增加ZMW孔的数目，从而来提升通量，ZMW孔的数目也从刚开始的3千个，一路增加到15万个、100万个和800万个，即将在2023年上市的Revio机型，单个芯片会包含2500万个ZMW孔。

当然这对于DNA分子较长的文库，和所需测序数据量较高的项目来说没有什么影响。但是实际中像细菌16S项目，全长的16S也只有1.5kb，或者转录组文库，反转录后的DNA分子长度平均也就1.5-2kb左右。就目前测序的酶读长，对于HiFi 15kb的文库数据，平均测序准确性可以达到QV30左右，即千分之一的错误率。如果文库较短，虽然可以获取更多的pass数，经过CCS矫正后具有较高的准确性，但是分子长度较短会使有效数据量变少，所以MAS-seq应运而生。

2. 文库构建过程

MAS-seq文库构建原理简单介绍如下：

构建样本的短片段文库，并将其平均分为多个独立的子文库，下图1是一个单细胞转录组文库构建过程，它分为了4个独立的子文库；

在每个文库中两端添加不同的接头，例如cDNA1文库的5'和3'分别添加AB接头，cDNA2文库的DNA分子5‘和3’分别添加B'C接头，cDNA3添加C'D接头，cDNA4添加D'E接头。再将这四个文库混成一个文库，由于BB'、CC'、DD'反向互补，连接生产一个较长片段文库；

目前PacBio官方提供8个子文库接头试剂盒，在15kb片段长度时就可以HiFi文库可以得到较高测序质量的数据，那么平均子文库的长度长度15/8约等于2kb；

图1. 单细胞转录组MAS-seq

3. MAS-seq数据拆分

PacBio测序仪可以支持Subreads类型下机数据和HiFi类型下机数据。两者之间的区别是，HiFi数据是将同一个ZMW孔中相同分子多次测到的Subreads经过合并纠错后生成的，具有较高的准确性，这一步转换称为CCS过程。在实际项目中，建议和服务商要求提供HiFi数据，因为这一转换过程非常消耗CPU资源。在CCS过程中，软件是不会去除MAS-seq数据中的接头序列的。为了将接头去除，获取实际DNA分子的碱基序列，可以使用软件Skera。拆分后的数据通常称为Segmented reads.

图2. Skera对MAS-seq数据拆分

5. 演示

官方提供了一组MAS-seq文库结构的HiFi数据，我们可以拿过来学习，了解具体的输入和输出数据格式，以及Skera软件的用法。

# download HiFi reads for MAS-Seq PBMCs run on Sequel IIe
wget https://downloads.pacbcloud.com/public/dataset/MAS-Seq/DATA-SQ2-PBMC_5kcells/0-CCS/m64476e_220618_014917.hifi_reads.bam

# download MAS adapter fasta
wget https://downloads.pacbcloud.com/public/dataset/MAS-Seq/REF-MAS_adapters/MAS-Seq_Adapter_v1/mas16_primers.fasta

# run skera split to generate segmented reads
skera split m64476e_220618_014917.hifi_reads.bam mas16_primers.fasta segmented.bam

其中mas16_primers.fasta为包含adapter序列文件，文件中接头序列的存放顺序，必须按照多个小片段文库连接生成最终长片段文库中，adapter从5‘到3’的顺序存放。

4. 输出文件

skera split 运行完成后会生成很多文件，它们包含不同信息。

文件名称	内容说明
segmented.bam	去除adapter后的测序数据
segmented.non_passing.bam	两端adapter不合理的序列
summary.csv	去除adapter后数据的统计信息
ligations.csv	5'和3'不同adapter组合的统计
read_lengths.csv	HiFi数据和Segmented reads长度统计信息

*summary.csv 输出文件格式说明:

Input Reads,2622891  #输入HiFI数据的总reads数吗
Segmented Reads (S-Reads),40131832  #被分割后得到segmented Reads总数
Mean Length of S-Reads,672 #segmented reads平均长度
Percentage of Reads with Full Array,86.3247  #拆分时一条HiFi reads得到的segmented reads数与文库设计一致的比例，本实例中包含15个segmented reads的HiFi reads数目占输入条数的比例
Mean Array Size (Concatenation Factor),15  #文库构建时完整HiFi reads理论包含segment数目

拆分数据时只输出按照adapter顺序连接的segment reads，例如官方提供的测试数据包含16个adapter，按照1-2，2‘-3，3’-4 ... 15‘-16 的方式加入到小片段文库中，然后通过2-2‘，3-3’ ... 15-15'粘性末端连接的方式将15条segment reads连接成一条长序列用于测序。由于在小片段文库中添加adapter是独立的，所以不应该产生1-3’、3-6‘这样连接方式，在拆分时也不回输出，下例中可以看出顺序连接的数目也是最多的。

sed 's/,/\t/g' m64476e_220618_014917.skera.ligations.csv|sort -k3,3rn | awk '{sum+=$3;print $0,sum}'|less -S

1       2       2579711 2579711
2       3       2569629 5149340
3       4       2555360 7704700
4       5       2548285 10252985
5       6       2535822 12788807
6       7       2522046 15310853
7       8       2509664 17820517
8       9       2502244 20322761
9       10      2495860 22818621
0       1       2494407 25313028
10      11      2490646 27803674
11      12      2486362 30290036
12      13      2481755 32771791
13      14      2478636 35250427
14      15      2475053 37725480
15      16      2406352 40131832
#只输出上面连接的segmented reads
2       4       8220 40140052
7       16      3157 40143209
6       16      2915 40146124
1       3       2885 40149009
5       16      2327 40151336
...
...

segment.bam文件中包含去除adapter的数据，Skera还会在每个序列后面添加很多tag以表示不同的信息。

BAM tag名	类型	含义	举例
di	i	segment的编号	di:i:0
qs	i	segment在原始HiFI数据中的起始位置	qs:i:16
qe	i	segment在原始HiFI数据中的终止位置	qe:i:450
dl	i	5'端adapter在adatper fasta文件中的索引	dl:i:0
dr	i	3'端adapter在adapter fasta文件中的索引	dr:i:1
ds	b	binary json，用于将segment还原为HiFi数据	ds:b:10,21,23

另外，被分割生成的segment reads可以反向生成连接到一起的HiFi数据，由于去掉了非顺序连接的segments reads，所以生成的undo bam会比原来的小：

skera undo *.skera.bam *undo.bam   #undo bam就是反向生成的HiFi数据的bam

MAS-seq技术原理（PacBio)