在SAM输出的结果中每一行都包括十二项通过Tab分隔(\t),缺失值使用’*'或者’0’代替,从左到右分别是:
1 QNAME,序列的名字(Read的名字)
2 FLAG, 概括出一个合适的标记,各个数字分别代表:
- 1 : 代表这个序列采用的是PE双端测序
- 2: 代表这个序列和参考序列完全匹配,没有插入缺失
- 4: 代表这个序列没有mapping到参考序列上
- 8: 这个序列是pair中的一个但是没有找到位点。代表这个序列的另一端序列没有比对到参考序列上,比如这条序列是R1,它对应的R2端序列没有比对到参考序列上
- 16:代表这个序列比对到参考序列的负链上,序列与参考序列反向互补
- 32 :代表这个序列对应的另一端序列比对到参考序列的负链上,与参考序列反向互补
- 64 : 代表这个序列是R1端序列, read1;
- 128 : 代表这个序列是R2端序列,read2;
- 256: 代表这个序列不是主要的比对,一条序列可能比对到参考序列的多个位置,只有一个是首要的比对位置,其他都是次要的
- 512: 代表这个序列在QC时失败了,被过滤不掉了(# 这个标签不常用)
- 1024: 代表这个序列是PCR重复序列(#这个标签不常用)
- 2048: 代表这个序列是补充的比对(#这个标签具体什么意思,没搞清楚,但是不常用)
如果flag值是0,那么说明测序为单端测序且这条read是primary line,一般是该read的最佳比对。这里说一下secondary alignment和supplementary alignment的区别,secondary alignment代表的是该条read比对到多个位置,该条记录为次优比对,在双端测序中,代表hardclip,而supplementary alignment标记为该条比对记录为该read的补充比对记录,关于补充比对记录其实就是嵌合比对,这个常见于三代测序结果中,因为三代测序的reads很长,在比对的时候,遇到大片段的缺失,基因融合或者转座(translocation)时,结果的bam文件中就会出现多条记录,这些记录是相互的补充,有可能只是这条read的一部分,我在二代测序的结果中很容易找到很多的hardclip记录,但是却没有找到supplementary alignment的记录,在三代的比对结果中就比较多了,以下是命令和一个截图:
假如说标记为以上列举出的数目,就可以直接推断出匹配的情况。假如说标记不是以上列举出的数字,比如说83=(64+16+2+1),就是这几种情况值和。
获取secendary alignment的条目
samtools view -f 256 xxx.bam | les
注意看cigar值那一列,就是红框的那一列,都是带有Hardclip的
#获取supplement alignment的条目
samtools view -f 2048 xxx.bam | les
在三代测序的比对结果中常见2048的supplementary alignment。使用samtools view -f/-F参数可以获取包含/不包含指定flag的内容的条目。
flag的64和128为在(Paired-end sequencing)双端测序中的read1和read2,而flag为16和32则代表双端测序中一条read反向互补,即代表比对到了负链上。
3 RNAME,参考序列的名字(染色体,contig)
4 POS,在参考序列上的位置(染色体/contig上的位置)
5 MAPQ, mapping qulity 越高则位点越独特
比对有时并不能完全确定一个短的序列来自参考序列的哪个位置,特别是对那些比较简单的序列。但是会给出一个值来显示这段序列来自某个位点的概率值,这个值就是mapping qulity。Mapping qulity的计算方法是:Q=-10log10p,Q是一个非负值,p是这个序列不来自这个位点的估计值。
假如说一条序列在某个参考序列上找到了两个位点,但是其中一个位点的Q明显大于另一个位点的Q值,这条序列来源于前一个位点的可能性就比较大。Q值的差距越大,这独特性越高。
6 CIGAR,代表比对结果的CIGAR字符串,如37M1D2M1I,这段字符的意思是37个匹配,1个参考序列上的删除,2个匹配,1个参考序列上的插入。M代表的是alignment match(可以是错配)
#standard cigar:
M match
I insertion
D deletion
#extended cigar
N gap
S substitution
H hard clipping
P padding
= sequence match
X sequence mismatch
Cigar值中的=为完全匹配,X为不完全匹配,M为最优比对,从biostar上看到的一种解释方法如下:
7 RNEXT, mate 序列所在参考序列的名称; 下一个片段比对上的参考序列的编号,没有另外的片段,这里是’*‘,同一个片段,用’=‘;
8 PNEXT, mate 序列在参考序列上的位置;下一个片段比对上的位置,如果不可用,此处为0;
9 TLEN,估计出的片段的长度,当mate 序列位于本序列上游时该值为负值。Template的长度,最左边得为正,最右边的为负,中间的不用定义正负,不分区段(single-segment)的比对上,或者不可用时,此处为0。另一个解释为:TLEN模板长度,如果reads对在同一条参考序列上,则模板长度为第4列减去第8列加上cigar值中的Match碱基数目。
10 SEQ,read的序列;序列片段的序列信息,如果不存储此类信息,此处为’*‘,注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度;
11 QUAL,ASCII码格式的序列质量,格式同FASTQ一样。这个指的是碱基质量值,它的算法如下:
这里p指得是碱基判断错误的概率,英文:base-calling error probabilities ,这里在加33之前的值叫Pred quality score
之所以要加33是为了对应ascii码中的可见字符,如果一个碱基的质量为0,那么对应的就是ascii码中的33,也就是!
12 可选的字段(field)
AS:i 匹配的得分
XS:i 第二好的匹配的得分
YS:i mate 序列匹配的得分
XN:i 在参考序列上模糊碱基的个数
XM:i 错配的个数
XO:i gap open的个数
XG:i gap 延伸的个数
NM:i 经过编辑的序列
YF:i 说明为什么这个序列被过滤的字符串
YT:Z
MD:Z 代表序列和参考序列错配的字符串
示例:
HWI-ST170:265:5:44:14178:183344#0 145 1 62421 37 63M1I35M 18 56843949 0 CCTGTATACATAGTAATCAAAGTGTACCACTGGTCGGTGTTTGTGTTCAGGCCCCTGTTGGGTAATGTGCATGTGAAGACCTCAGGTGGTATAGTTTTG CEE?@F@BE@GGEGFBHHEDEEEDEEBEDHHBGHGGFHHDFHHHGGGGFFFEEEHFHFGFHHHHHFHHHFHHHHGHGHEHHHHHHHHHFHHHHHHHHHH RG:Z:DU23M01_Duroc XT:A:U NM:i:4 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:3 XO:i:1 XG:i:1 MD:Z:20T22C1A52
HWI-ST170:264:5:61:3024:21492#0 113 1 62421 37 63M1I29M = 6885283 6822868 CCTGTATACATAGTAATCAAAGTGTACCACTGGTCGGTGTTTGTGTTCAGGCCCCTGTTGGGTAATGTGCATGTGAAGACCTCAGGTGGTATA @:;9AFGCHFHHHEGGGHDCADA?E@EEDAHFHHFFHHHFFHHHHHHHHHFHFHDHHHHHHHHGHGHHGHFHHHHHHHFHHHHHHHHHHHHHH RG:Z:DU23M01_Duroc XT:A:U NM:i:4 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:3 XO:i:1 XG:i:1 MD:Z:20T22C1A46
FCC1L2FACXX:3:2106:15923:93264 99 1 1073 0 100M = 1461 488 TGTGAAGGCCCCCTGCTCTGACTGTGTTAGTGTCCATTTCTCCTTTTACGGTTGTAGCAGTTGCCTTCTACATTGCGGGGATCCTGTATTGGGTGCATGT ___eceeegfggggdgiiifghii[degfhfgfdffhhhfhfghiiighiiH^`Vbgfffihhiiiihhddbdgfgccca][^bbbbbccbca[X^Y_b_ MD:Z:98A1 PG:Z:MarkDuplicatesRG:Z:SRR949625 NM:i:1 AS:i:98 XS:i:98
FCC1L2FACXX:3:2111:13731:89147 163 1 1073 7 100M = 1438 465 TGTGAAGGCCCCGTGCTCTGACTGTGTTAGTGTCCATTTCTCCTTTTACGGTTGTAGCAGTTGCCTTCTACATTGCGGGGATCCTGTATTGGGTGCATAT bbaeeeeefggggefhiiiihiiicgghhigdehhhiiihiieffhiihihbggdghihfgffhihihf`geed_cecac]accb]_bcccbc^a_bcbb MD:Z:12C87 PG:Z:MarkDuplicatesRG:Z:SRR949625 NM:i:1 AS:i:95 XS:i:95
一些重要的概念
首先是关于模板和read的概念
Template: A DNA/RNA sequence part of which is sequenced on a sequencing machine or assembled from
raw sequences.
Segment: A contiguous sequence or subsequence.
Read: A raw sequence that comes off a sequencing machine. A read may consist of multiple segments. For
sequencing data, reads are indexed by the order in which they are sequenced
然后是线性比对(linear alignment)和嵌合比对(chimeric alignment)的概念
Linear alignment: An alignment of a read to a single reference sequence that may include insertions,
deletions, skips and clipping, but may not include direction changes (i.e., one portion of the alignment
on forward strand and another portion of alignment on reverse strand). A linear alignment can be
represented in a single SAM record.
Chimeric alignment: An alignment of a read that cannot be represented as a linear alignment. A chimeric
alignment is represented as a set of linear alignments that do not have large overlaps. Typically, one
of the linear alignments in a chimeric alignment is considered the “representative” alignment, and the
others are called “supplementary” and are distinguished by the supplementary alignment flag. All the
SAM records in a chimeric alignment have the same QNAME and the same values for 0x40 and 0x80
flags (see Section 1.4). The decision regarding which linear alignment is representative is arbitrary.
Read alignment: A linear alignment or a chimeric alignment that is the complete representation of the
alignment of the read.
Multiple mapping: The correct placement of a read may be ambiguous, e.g., due to repeats. In this case,
there may be multiple read alignments for the same read. One of these alignments is considered
primary. All the other alignments have the secondary alignment flag set in the SAM records that
represent them. All the SAM records have the same QNAME and the same values for 0x40 and 0x80
flags. Typically the alignment designated primary is the best alignment, but the decision may be
arbitrary.3
线性比对和嵌合比对:实际上,一条read的比对情况有时候会很复杂,例如一个read可能比对到多个位置,因为基因组非常庞大,存在很多的repeat区域,而二代测序的read一般比较短,在200bp以内,所以这种时候,可能在bam文件就会出现多条记录该read比对情况的记录,它们拥有相同的id,其中的一个较好的比对被认为是"典型的"比对结果,而其他的都被称为“补充的”比对结果。而且就算比对到同一个位置,还有可能有方向性问题,这是也会有两条记录存在,而没有方向问题,单一的最优比对的read就被叫做线性比对,具体表现在bam文件中将只有一条关于该read的比对记录(注意是read而不是paired reads)。
Phred scale
Phred scale: Given a probability 0 < p ≤ 1, the phred scale of p equals −10 log10 p, rounded to the closest
integer.
这个值是在二代测序文件中常常出现的Phred 值,这个值是对p值的一个log变换,可以更加直观的反应p值,Phred值越大,代表p值越小,例如Phred值为10,那么p值则为0.1;Phred为20时,p值为0.01 。
1-base和0-base文件
根据samtools的说明,1base的文件主要有 SAM,VCF,GFF,Wiggle文件,而0base的文件主要是BAM,BCF,BED这样的文件。
1-based coordinate system: A coordinate system where the first base of a sequence is one. In this coordinate
system, a region is specified by a closed interval. For example, the region between the 3rd
and the 7th bases inclusive is [3, 7]. The SAM, VCF, GFF and Wiggle formats are using the 1-based
coordinate system.
0-based coordinate system: A coordinate system where the first base of a sequence is zero. In this
coordinate system, a region is specified by a half-closed-half-open interval. For example, the region
between the 3rd and the 7th bases inclusive is [2, 7). The BAM, BCFv2, BED, and PSL formats are
using the 0-based coordinate system.
bam文件的header
找到一个bam文件并输入如下命令查看title信息
samtools view *.bam -H|less
可以看到bam文件的header有如下几个部分组成
- @HD 标题行(The header line)主要用于版本号,是否排序等标注
- @SQ 参考序列字典(Reference sequence dictionary)注意它们的顺序决定了排序的优先级
- @RQ read group主要记录了一文件的描述信息,例如样本信息,文库编号等
- @PG program主要记录比对的软件等代码信息
- @CO comments一行评论
//www.greatytc.com/p/c48c36affff7
https://samtools.github.io/hts-specs/SAMv1.pdf
https://genome.sph.umich.edu/wiki/SAM
https://en.wikipedia.org/wiki/SAM_(file_format)
https://www.biostars.org/p/60765/
https://www.cnblogs.com/emanlee/p/4316581.html
https://blog.csdn.net/qq_35696312/article/details/101760397