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}
注意:
- 这种未比对的bam必须是按第一列排序的。所有来自同一个ZMW hole的subreads放在一起,且按qStart排序,每个ZMW再按数字排序。特别是从subreads bam文件中随机取subreads分析时一定要用samtools sort -n处理一下
- 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的转化如下:
Frames | Encoding |
0, 1, ... 63 | 0, 1, ... 63 |
64, 66, ... 190 | 64, 65, ... 127 |
192, 196, ... 444 | 128, 129, ... 191 |
448, 456, ... 952 | 192, 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