Pacbio Sequel两种bam文件解析

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

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

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

@HD     VN:1.5  SO:unknown      pb: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? 代表序列和参考序列错配的字符串

附加参考材料:

[1] https://samtools.github.io/hts-specs/SAMv1.pdf

[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

 

发表评论

匿名网友