热点新闻
生信分析常见文件——BAM文件
2024-10-15 14:37  浏览:4062  搜索引擎搜索“手机全球会展网”
温馨提示:信息一旦丢失不一定找得到,请务必收藏信息以备急用!本站所有信息均是注册会员发布如遇到侵权请联系文章中的联系方式或客服删除!
联系我时,请说明是在手机全球会展网看到的信息,谢谢。
展会发布 展会网站大全 报名观展合作 软文发布

比对软件将质控后的fq格式文件与参考基因组进行比对后得到SAM文件,首先由于sam文件很大,mapping结果乱序不可直接用,将sam文件转为bam格式,再对bam文件进行排序sorting

一个非常好用的工具包——Samtools。

SAM(The Sequence Alignment / Map format)文件是BWA比对软件输出的纯文本文件,BAM(B取自binary)文件是SAM文件的(压缩为1/6)。BAM文件由两部分组成,头部区(header)和主体区(record),都以tab分列。

一、头部区简要介绍

头部区:以’@'开始,体现了比对的一些总体信息。比如比对的SAM格式版本,比对的参考序列,比对使用的软件等。

@HD VN:1.0 SO:unsorted (排序类型)

头部区第一行:VN是格式版本;SO表示比对排序的类型,有unknown(default),unsorted,queryname和coordinate几种。samtools软件在进行行排序后不能自动更新bam文件的SO值,而picard却可以。

@SQ SN:contig1 LN:9401 (序列ID及长度)

参考序列名,这些参考序列决定了比对结果sort的顺序,SN是参考序列名;LN是参考序列长度;每个参考序列为一行。

例如:@SQ SN:NC_000067.6 LN:195471971

@RG ID(read group):sample01 (样品基本信息)

1个sample的测序结果为1个Read Group;该sample可以有多个library的测序结果,可以利用bwa mem -R 加上去这些信息。

例如:@RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq

ID:样品的ID号 SM:样品名 LB:文库名 PU:测序以 PL:测序平台

这些信息可以在形成sam文件时加入,ID是必须要有的后面是否添加看分析要求

@PG ID:bowtie2 PN:bowtie2 VN:2.0.0-beta7 (比对所使用的软件及版本)

例如:@PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa sampe -a 400 -f ZX1.sam -r @RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq …/0_Reference/Reference_Sequence.fa ZX_HQ_clean_R1.fq.sai ZX_HQ_clean_R2.fq.sai …/2_HQData/ZX_HQ_clean_R1.fq …/2_HQData/ZX_HQ_clean_R2.fq

这里的ID是bwa,PN是bwa,VN是0.7.12-r1039版本。CL可以认为是运行程序@RG是上面RG表示的内容,后面是程序内容,这里的@GR内容是可以自己在运行程序是加入的。

二、主体部分介绍

主体区:比对结果,每一个比对结果是一行,有11个主列和1个可选列。

(如果某一列为“0”或“*”表示这一列没有信息)













(1)第一列:QNAME:比对的序列名称(reads ID)。

进行reads比对时通常表示reads的名字,如果这条reads比对到多条序列或比对到这条序列的多个位置,相同名字会出现多次。如果是pair-end reads,相同名字会出现2次,分别表示来自于R1文件的reads和R2文件的reads,如果其matepair reads也比对2个位置,也会出现2次,则相同名字共出现4次,如果一条reads也比对2个位置,则其matepair比对1个位置,则共出现3次,如果其matepair reads没有比对上序列也会出现1次(第三列显示“*”),所以pair-end测序,R1文件和R2文件同时mapping,相同reads的id最少出现2次。

(2)第二列:FLAG:比对信息位

数值结果如下:

1(1)该read是成对的paired reads中的一个

2(10)paired reads中每个都正确比对到参考序列上

4(100)该read没比对到参考序列上

8(1000)与该read成对的matepair read没有比对到参考序列上

16(10000)该read其反向互补序列能够比对到参考序列

32(100000)与该read成对的matepair read其反向互补序列能够比对到参考序列

64(1000000)在paired reads中,该read是与参考序列比对的第一条

128(10000000)在paired reads中,该read是与参考序列比对的第二条

256(100000000)该read是次优的比对结果

512(1000000000)该read没有通过质量控制

1024(10000000000)由于PCR或测序错误产生的重复reads

2048(100000000000)补充匹配的read

具体的flag值的解释,可以参考samtools软件提供的结果

其中的samtools flags用法可提供flag值的查找结果

Usage: samtools flags INT|STR[,…]

如果某一个数值不是下面的任意值,那么那个数值就是下面这些数里面几个的和。

例如:samtools flags 10→0xa 10 PROPER_PAIR,MUNMAP(10=2+8)

          samtools flags 12→0xc 12 UNMAP,MUNMAP(12=4+8)

(3)第三列:RNAME:染色体

表示read比对的那条序列的序列名称(名称与头部的@SQ相对应),如果这列是“*”,可以认为这条read没有比对上的序列,则这一行的第四,五,八,九 列是“0”,第六,七列与该列是相同的表示方法。

(4)第四列:POS:表示比对的起始位置,以1开始计数,如果没有比对上则显示为0。

(5)第五列:MAPQ:表示为mapping的质量值,mapping Quality。

比如说Q30(错配率为0.001),Q20(错配率为0.01),计算公式为:-10logP{错比概率} 。如果值为255表示mapping值是不可用的,如果是unmapped read则MAPQ为0,一般在使用bwa mem生成的sam文件,第五列为60表示mapping率最高,一般结果是这一列的数值是从0到60,且0和60这两个数字出现次数最多。

(6)第六列:CIGAR:reads mapping到第三列序列的mapping状态,

对于mapping状态可分为以下几类:

M(匹配):alignment match (can be a sequence match or mismatch)

表示read可mapping到第三列的序列上,则read的碱基序列与第三列的序列碱基相同,表示正常的mapping结果,M表示完全匹配,但是无论reads与序列的正确匹配或是错误匹配该位置都显示为M

I(插入):insertion to the reference

表示read的碱基序列相对于第三列的RNAME序列,有碱基的插入

D(删除):deletion from the reference

表示read的碱基序列相对于第三列的RNAME序列,有碱基的删除

N:skipped region from the reference    表示可变剪接位置

P:padding (silent deletion from padded reference)

S:soft clipping (clipped sequences present in SEQ)  这部分没比对上但保留了

H:hard clipping (clipped sequences NOT present in SEQ) 这部分没比对上且未保留

clipped均表示一条read的序列被分开,之所以被分开,是因为read的一部分序列能匹配到第三列的RNAME序列上,而被分开的那部分不能匹配到RNAME序列上。

"="表示正确匹配到序列上

"X"表示错误匹配到序列上

而H只出现在一条read的前端或末端,但不会出现在中间,S一般会和H成对出现,当有H出现时,一定会有一个与之对应的S出现

比如3S6M1P1I4M,前三个碱基被剪切去除了,然后6个比对上了,然后打开了一个缺口,有一个碱基插入,最后是4个比对上了,是按照顺序的;







(7)第七列:MRNM

这条reads第二次比对的位置,=表示参考序列与reads一模一样,*表示没有完全一模一样的参考序列。

(8)第八列:MPOS该列表示与该reads对应的mate pair reads的比对位置(即mate),若无mate,则为0。

(9)第九列:ISIZE序列模板长度,如果同一个片段都比对上了同一个参考序列,为最左边的碱基位置到最右边的碱基位置(左为正,右为负)。当mate 序列位于本序列上游时该值为负值。不可用时,为0。

(10)第十列:read的序列。

(11)第十一列:ASCII码格式的序列质量。格式同FASTQ一样。其中1、10、11合起来就是fq格式文件。

(12)第十二列:可选的区域。格式类似AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:35T0 YT:Z:UU

处理bam文件的主要生信软件有

bwa,bowtie2,samtools,bedtools等

可以看mapping等多方面结果和统计,bedtools工具中genomeCoverageBed的功能是:Compute the coverage of a feature file among a genome

三、查看BAM文件

BAM文件是二进制格式,没办法通过文本的形式直接打开,要用samtools的view功能在终端上进行查看。

$ samtools view in.bam  #重头开始查看(不建议)

如果不想从头开始看,希望快速地跳转到基因组的其它位置上,比如chr22染色体,那么可以先用samtools index生成BAM文件的索引(如果已经有索引文件则不需该步骤),然后进行如下操作:

$ samtools index in.bam  # 生成in.bam的索引文件in.bam.bai

$ samtools view in.bam chr22            # 跳转到chr22染色体

$ samtools view in.bam chr22:16050103  # 跳转到chr22:16050103位置

$ samtools view in.bam chr22:16050103-16050103  # 只查看该位置







参考:bam格式文件详解-CSDN博客

           bam文件的理解 - 简书 (jianshu.com)

           第5节-理解并操作BAM文件 · 从零开始完整学习全基因组测序数据分析 (gitbooks.io)

发布人:1570****    IP:124.223.189***     举报/删稿
展会推荐
让朕来说2句
评论
收藏
点赞
转发