本文链接:https://blog.csdn.net/genome_denovo/article/details/78712972
第二章:生物信息分析
第一节:解读sam格式文件
1,SAM文件格式介绍
SAM(The Sequence Alignment / Map format)格式,即序列比对文件的格式,详细介绍文档:http://samtools.github.io/hts-specs/SAMv1.pdf
SAM文件由两部分组成,头部区和主体区,都以tab分列。
头部区:以’@'开始,体现了比对的一些总体信息。比如比对的SAM格式版本,比对的参考序列,比对使用的软件等。
主体区:比对结果,每一个比对结果是一行,有11个主列和一个可选列。
2,头部区简要介绍
@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:sample01 (样品基本信息)
Read Group。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内容是可以自己在运行程序是加入的
3,主体部分介绍
主体部分有11个主列和1个可选列
QNAME 比对的序列名称 例如:M04650:84:000000000-B837R:1:1101:22699:1759(一条测序reads的名称)
FLAG Bwise FLAG(表明比对类型:paring,strand,mate strand等) 例如:99
RENAME 比对上的参考序列名 例如:NC_000075.6
POS 1-Based的比对上的最左边的定位 例如:124057649
MAPQ 比对质量 例如:60
CIGAR Extended CIGAR string(操作符:MIDNSHP)比对结果信息;匹配碱基数,可变剪接等 例如:87M
MRNM 相匹配的另外一条序列,比对上的参考序列名 例如:=
MPOS 1-Based leftmost Mate Position (相比于MRNM列来讲意思和POS差不多) 例如:124057667
ISIZE 插入片段长度 例如:200
SEQ 和参考序列在同一个链上比对的序列(若比对结果在负义链上,则序列是其反向重复序列,反向互补序列) 例如:ATTACTTGGCTGCT
QUAL 比对序列的质量(ASCII-33=Phred base quality)reads碱基质量值 例如:-8CCCGFCCCF7@E-
可选的列 以TAG:TYPE:VALUE的形式提供额外的信息
4,对于每一列内容的详细注解
(如果某一列为“0”或“*”表示这一列没有信息)
第一列:QNAME
进行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次。
第二列: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(Version: 1.3.1)
其中的samtools flags用法可提供flag值的查找结果
About: Convert between textual and numeric flag representation
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)
具体的flag值的解释,也可参考如下网站:https://broadinstitute.github.io/picard/explain-flags.html
或者在必应当中搜索flag sam点击Explain SAM Flags-GitHub Pages进入该网页,也可以输入组合flag数值会出现所存在的意思
第三列:RNAME
表示read比对的那条序列的序列名称(名称与头部的@SQ相对应),如果这列是“*”,可以认为这条read没有比对上的序列,则这一行的第四,五,八,九 列是“0”,第六,七列与该列是相同的表示方法
第四列:POS
表示read比对到RNAME这条序列的最左边的位置,如果该read能够完全比对到这条序列(CIGAR string为M)则这个位置是read的第一个碱基比对的位置,如果该read的反向互补序列比对到这条序列,则这个位置是read的反向互补序列的第一个碱基比对的位置,所以无论该read是正向比对到该序列,或是其反向互补序列比对到该序列,比对结果均是最左端的比对位置
第五列:MAPQ
表示为mapping的质量值,mapping Quality, It equals -10log10Pr{mapping position is wrong}, rounded to the nearest integer, A value 255 indicates that the mapping quality is not available. 该值的计算方法是mapping的错误率的-10log10值,之后四舍五入得到的整数,如果值为255表示mapping值是不可用的,如果是unmapped read则MAPQ为0,一般在使用bwa mem或bwa aln(bwa 0.7.12-r1039版本)生成的sam文件,第五列为60表示mapping率最高,一般结果是这一列的数值是从0到60,且0和60这两个数字出现次数最多
第六列:CIGAR
CIGAR string,可以理解为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出现
例如:
162M89S
162H89M
149M102S
149H102M
40S211M
20M1D20M211H
S可以单独出现,而H必须有与之对应的S出现时才可能出现,不可在相同第一列的情况下单独出现
N:如果是mRNA-to-genome,N出现的位置代表内含子,其它比对形式出现N时则没有具体解释
M/I/S/=/X:这些数值的加和等于第10列SEQ的长度
第七列:MRNM
这条reads第二次比对的位置,在利用bwa mem产生sam文件时,如果该列是“”而
第3列RNAME不是“”则表示该reads比对到第3列显示序列名的序列上,而没有比对到其他位置,在利用bwa aln及bwa sampe比对生成的sam文件,如果和上述情况相同,则第7列为“=”,上述情况均表示该reads只比对到这一个位置
如果第3列RNAME和第7列MRNM都为“*”,则说明这条reads没有匹配上的序列,如果这条reads匹配两个序列,则第一个序列的名称出现在第3列,而第二个序列的名称出现在第7列
第八列:MPOS
该列表示与该reads对应的mate pair reads的比对位置,如果这对pair-end reads比对到同一条reference序列上,在sam文件中reads的id出现2次,Read1比对的第4列等于Read2比对的第8列。同样Read1比对的第8列等于Read2比对的第4列。例如:
第1列(Read id)····第4列(Read1比对位置)····第8列(mate-pair reads比对位置)
22699:1759····124057649····124057667
22699:1759····124057667····124057649
相同的reads id一个来自Read1文件,一个来自Read2文件,第4列和第8列是对应的
第九列:ISIZE
TLEN:signed observed Template LENgth (可以理解为文库插入片段长度)
如果R1端的read和R2端的read能够mapping到同一条Reference序列上(即第三列RNAME相同),则该列的值表示第8列减去第4列加上第6列的值,R1端和R2端相同id的reads其第九列值相同,但该值为一正一负,R1文件的reads和R2文件的reads,相同id的reads要相对来看。在进行该第列值的计算时,如果取第6列的数值,一定要取出现M的值,S或H的值不能取。
the unisgned observed template length equals the number of base from the leftmost mapped base to the rightmost mappedbase. Theleftmost segment has a plus sign and the rightmost has a minus sign
处理bam文件的主要生信软件有
bwa,bowtie2,samtools,bedtools等
可以看mapping等多方面结果和统计,bedtools工具中genomeCoverageBed的功能是:Compute the coverage of a feature file among a genome
————————————————
版权声明:本文为CSDN博主「genome_denovo」的原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接及本声明。
原文链接:https://blog.csdn.net/genome_denovo/article/details/78712972