sam格式很精炼，几乎包含了比对的所有信息，我们平常用到的信息很少，但特殊情况下，我们会用到一些较为生僻的信息，关于这些信息sam官方文档的介绍比较精简，直接看估计很难看懂。
今天要介绍的是如何通过bam文件统计比对的indel和mismatch信息
首先要介绍一个非常重要的概念--编辑距离
定义：从字符串a变到字符串b，所需要的最少的操作步骤（插入，删除，更改）为两个字符串之间的编辑距离。
这也是sam文档中对NM这个tag的定义。
编辑距离是对两个字符串相似度的度量（参见文章：Edit Distance）
举个栗子：两个字符串“eeba”和“abca”的编辑距离是多少？
根据定义，通过三个步骤：1.将e变为a 2.删除e 3.添加c，我们可以将“eeba”变为“abca”
所以，“eeba”和“abca”之间的编辑距离为3
回到序列比对的问题上
下面是常见的二代比对到ref的结果（bwa）：
D00360:96:H2YLYBCXX:1:2110:18364:84053 353 seq1_len154_cov5 1 1 92S59M8I17M1D6M1D67M seq30532_len2134_cov76 1 0 AAAAAAAAAAAAAAAAAAAAAAAACCCTGTCTCTAATAAAATACAAACAATTAGCCGGGCATGGTGGCACGCGCCTTTAGTCCCAGCTACTAGGGAGGCTGAGGCAGGGGAATTGTTTGAACCCGGGAGGTGGAGGTTGCAGTGAGCGGAGTTTTTTTCACTGCACTCCAGCCTGGTGACAAATCAAAAATCCATTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAACAACAAA
DDDDDHIIIIIIII<EHHII?EHH0001111111<11<DEH1D1<FH1D<<<C<@GEHD</<11<101<1D<<C<E0D11<<1<D?1F1CC1DE110C0D1011100////0DD<[email protected]=FGHCDHH<FG<D0<<<EF?CE<00<<0<D//0;<:D/////;///////;8F.;/.8.8......88.9........-8BBGADHIIHD?>D?HH<,>=HHDD,5CHDCDHD><,,,--8----8-8-- NM:i:25 MD:Z:16A21C16C0A3T15^A6^G1G0T0G3C2T0G1C41A4A2A3 AS:i:49 XS:i:42 SA:Z:seq13646_len513_cov63,125,-,103S125M21S,1,11; RG:Z:chr22
这是ref序列
>seq1_len154_cov5 GGGAGGCTGAGGCAGGAGAATTGTTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCAAGATTGCACTCCAGCCTGGATGACAAGAGTGAAACTCTGTCTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
cigar字段包含了序列比对的简化信息，M（匹配比对，包含match和mismatch），=（纯match），X（纯mismatch），I（插入），D（删除），还有N、P和S、H。（注：目前只在blasr比对结果中见过=和X）
根据cigar字段可以统计indel信息，但是无法统计mismatch。
这个时候就可以用到NM tag了，mismatch = NM – I - D = 25 – 8 – 1 – 1 = 15
有兴趣的可以对着cigar数一遍，下面是我无聊数着的结果，也是15个
使用python的pysam模块可以很容易的提取出每个比对结果的NM信息
参见之前的帖子：pysam - 多种数据格式读写与处理模块(python)
时间： 2024-10-16 03:59:40