pacbio目前有两种主流的测序平台，RSII和Sequel，后者是前者的升级版。

pacbio sequel下机是bam格式的reads文件，它和reads比对到参考基因组上生成的bam文件，内容有差异，但格式一致。

1. pacbio sequel下机的bam格式reads文件header部分如下：

@HD     VN:1.5  SO:unknownpb:3.0.3

@RG     ID:ceba59d4     PL:PACBIO       DS:READTYPE=SUBREAD;Ipd:CodecV1=ip;PulseWidth:CodecV1=pw;BINDINGKIT=100-862-200;SEQUENCINGKIT=100-861-800;BASECALLERVERSION=5.0.0.6236;FRAMERATEHZ=80.000000    PU:m54191_170909_212150 PM:SEQUEL

@PG     ID:baz2bam      PN:baz2bam      VN:5.0.0.6236   CL:/opt/pacbio/ppa-5.0.0/bin/baz2bam/data/pa/m54191_170909_212150.baz -o /data/pa/m54191_170909_212150 --metadata/data/pa/.m54191_170909_212150.metadata.xml -j 12 -b 12 --progress --silent--minSubLength 50 --minSnr 3.750000 --adapters /data/pa/m54191_170909_212150.adapters.fasta

@PG     ID:bazFormat    PN:bazformat    VN:1.3.0

@PG     ID:bazwriter    PN:bazwriter    VN:5.0.0.6236

# @RG中的DS是describe的意思。BINDINGKIT 和SEQUENCINGKIT 在Pacbio RSII和Pacbio Sequel中不同。PU即platform unit, 是Pacbio moive名，一次moive会用到多个零模波导孔(ZMW)，每个ZMW通常会有多条Subreads，多条subreads对应1个CCS(如果测序片段短于ZMW一次能测到的长读，则会产生哑铃型结构，循环来测这一段)。

2. body部分如下：

第一列Qname如下：

m54191_170909_212150/4391389/0_5019

#命名规则如下：

{movieName}/{holeNumber}/{qStart}_{qEnd}

注意：

1.这种未比对的bam必须是按第一列排序的。所有来自同一个ZMW hole的subreads放在一起，且按qStart排序，每个ZMW再按数字排序。特别是从subreads bam文件中随机取subreads分析时一定要用samtools  sort -n处理一下

2. qStart是从0开始的，qEnd – qStart =第10列的碱基序列长度。

第二列是flag值：

如果是Subreads则全部为4表示read unmapped。在三代长reads比对结果中会出现两种：0或16

第三列是ref值：

如果是Subreads则全部为*

第四列是reads在参考序列上的位置

如果是Subreads则全部为0

第五列比对质量值

如果是Subreads则全部为255。

Mapping qulity的计算方法是：Q=-10log10p，Q是一个非负值，p是这个序列来自这个位点的估计值。

假如说一条序列在某个参考序列上找到了两个位点，但是其中一个位点的Q明显大于另一个位点的Q值，这条序列来源于前一个位点的可能性就比较大。Q值的差距越大，这独特性越高[摘自http://www.bbioo.com/lifesciences/40-113305-1.html]。

第六列是CIGAR字符串

如果是Subreads则全部为*。如果是三代reads比对的结果与二代不同，不再用M字符(在二代的bam中它可以容许错配)，三代中采用更为精确的X(它表示与ref不同)和=(它表示与ref相同)，如果CIGAR中检测到M字符Pacbio software将会中止运行。D表示deletion，I表示insertion。整个CIGAR字符串首尾的H和S表示硬切和软切(soft-clippedbase)

第七列是mate序列比对到的参考序列的名称

因为三代没有PE，Subreads bam和aligned bam中全部为*

第八列是mate序列比对到参考序列上的位置

同第七列，Subreads  bam和aligned  bam全部为0

第九列为估计出的插入片段长度，当mate序列位于本序列上游时该值为负值。

同第七列，Subreads  bam和aligned  bam全部为0

第十列为read序列

第十一列是ASCII码格式的序列质量值

如果是Subreads则每个碱基质量值为！

以下各列格式如下

列名简写(两个小写字母)：数据类型：数据值

#数据类型：A:可打印字符，Z:可打印字符串，H: Hex字符，f:单精度浮点数[http://blog.csdn.net/arnoyshell/article/details/9383555]

第十二列是cx列，表示subread local context

LocalContextFlags对应情况如下：

ADAPTER_BEFORE = 1,

ADAPTER_AFTER  = 2,

BARCODE_BEFORE = 4,

BARCODE_AFTER  = 8,

FORWARD_PASS   = 16,

REVERSE_PASS   = 32

因为pacbio sequel下机数据已经切除了两边接头，可以不用关心这列

第十三列为ip列使用的是IPD (raw frames or codec V1)转化过的质量值。IPD的转化如下：FramesEncoding

0, 1, ... 630, 1, ... 63

64, 66, ... 19064, 65, ... 127

192, 196, ... 444128, 129, ... 191

448, 456, ... 952192, 193, ... 255

即第一个0-63直接编码，63-2*63每两个值合一，3*63-4*63每三个值合一[http://pacbiofileformats.readthedocs.io/en/3.0/BAM.html]

第十四列为np列表示NumPasses

第十五列为pw列表示PulseWidth

第十六列为qe列，表示query end

第十七列为qs列，表示query start

第十八列为rq列，用[0-1]编码的期望精确度

第十九列为sn列，4个float数依次表示ACGT平均信噪比

第二十列为zm列，表示零模波导孔号

第二十一列为ReadGroup列

另外可选的有

AS:i?匹配的得分

XS:i?第二好的匹配的得分

YS:i? mate序列匹配的得分

XN:i?在参考序列上模糊碱基的个数

XM:i?错配的个数

XO:i? gap open的个数

XG:i? gap延伸的个数

NM:i?经过编辑的序列，read string转换成reference string需要的最少核苷酸:插入/缺失/替换

YF:i?说明为什么这个序列被过滤的字符串

MD:Z?代表序列和参考序列错配的字符串

附加参考材料：

[2] http://pacbiofileformats.readthedocs.io/en/3.0/BAM.html

[3] https://github.com/PacificBiosciences/PacBioFileFormats/blob/3.0/BAM.rst

[4] http://albiorix.bioenv.gu.se:8081/smrtanalysis/doc/bioinformatics-tools/pbsamtools/doc/index.html

转载本文请联系原作者获取授权，同时请注明本文来自卢锐科学网博客。

链接地址：http://blog.sciencenet.cn/blog-2970729-1086176.html

上一篇：Speedseq的安装和使用

下一篇：比对算法的原理及代码实现