生物信息数据格式：sam,bam格式

最新推荐文章于 2024-06-23 09:48:02 发布

sunchengquan

最新推荐文章于 2024-06-23 09:48:02 发布

阅读量5.1k

点赞数 3

分类专栏： bioinformation

本文链接：https://blog.csdn.net/sunchengquan/article/details/85095355

版权

bioinformation 专栏收录该内容

21 篇文章 52 订阅

订阅专栏

文章目录

数据获取
格式说明
实例演练

数据获取

首先安装bowtie短序列比对软件

bam和sam文件主要是bwa、bowtie、tophat等序列比对工具产生的

wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip

unzip bowtie2-2.3.4.3-linux-x86_64.zip
ln -s ~/local/app/bowtie2-2.3.4.3-linux-x86_64/bowtie2 ~/bin

生成sam,bam文件

cd ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads/

bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq 2>alignment.log >tmp.sam


bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq 2>alignment.log |samtools view -bS  >tmp.bam

bam 是sam二进制格式

格式说明

sam(Sequence Alignment Mapping) 序列比对映射

sam分为两部分，注释信息（header section）和比对结果部分（alignment section）

注释信息可有可无，都是以@开头，用不同的tag表示不同的信息，主要有

@HD，说明符合标准的版本、对比序列的排列顺序；
@SQ，参考序列说明；例如 @SQ SN:chrY LN:57227415
@RG，比对上的序列（read）说明；
@PG，使用的程序说明；例如 @PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa mem /home/sunchengquan/data/database/hg38.fa IonXpress_001.fq
@CO，任意的说明信息。

比对结果部分（alignment section），每一行表示一个片段（segment）的比对信息，包括11个必须的字段（mandatory fields）和一个可选的字段，字段之间用tag分割。必须的字段有11个，顺序固定，不可用时，根据字段定义，可以为’0‘或者’*‘

QNAME,序列的名字（Read的名字）
FLAG, 概括出一个合适的标记，各个数字分别代表
1 read有多种测序数据，一般理解为有双端测序数据 read paired
2 比对结果是一个pair-end比对，一正一负完美的比对上，read mapped in proper pair
4 这条read没有比对上，read unmapped
8 这个序列是pair中的一个但是没有比对上 mate unmapped
16 序列与参考序列反向互补 read reverse strand
32 这个序列在pair-end中的的mate序列与参考序列反响互补 mate reverse strand
64 序列是 mate 1 first in pair
128 序列是 mate 2 second in pair
256 第二次比对 not primary alignment
…
RNAME,参考序列的名字（染色体）
POS,在参考序列上的位置（染色体上的位置）
MAPQ, mapping qulity 越高则位点越独特
bowtie2有时并不能完全确定一个短的序列来自参考序列的哪个位置，特别是对那些比较简单的序列。但是bowtie2会给出一个值来显示这个段序列来自某个位点的概率值，这个值就是mapping qulity。Mapping qulity的计算方法是：Q=-10log10p，Q是一个非负值，p是这个序列不来自这个位点的估计值。
假如说一条序列在某个参考序列上找到了两个位点，但是其中一个位点的Q明显大于另一个位点的Q值，这条序列来源于前一个位点的可能性就比较大。Q值的差距越大，这独特性越高。
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
RNEXT, mate 序列所在参考序列的名称; 下一个片段比对上的参考序列的编号，没有另外的片段，这里是’*‘，同一个片段，用’=‘；
PNEXT, mate 序列在参考序列上的位置;下一个片段比对上的位置，如果不可用，此处为0；
TLEN，估计出的插入片段的长度，当mate 序列位于本序列上游时该值为负值。Template的长度，最左边得为正，最右边的为负，中间的不用定义正负，不分区段（single-segment)的比对上，或者不可用时，此处为0；
SEQ，read的序列;序列片段的序列信息，如果不存储此类信息，此处为’*‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度；
QUAL，ASCII码格式的序列质量;序列的质量信息，格式同FASTQ一样。
可选的字段（field）
NM:i 经过编辑的序列
MD:Z 代表序列和参考序列错配的字符串
AS:i 匹配的得分

insert size示意图：

在这里插入图片描述

insertsize = mate位置B1（第8列）+mate长度（B2-B1）- 比对起始位置A1（4列）

FLAG含义查询网站： http://broadinstitute.github.io/picard/explain-flags.html ，假如说标记为以上列举出的数目，就可以直接推断出匹配的情况。假如说标记不是以上列举出的数字，比如说 83=（64+16+2+1），就是这几种情况值和。

实例演练

统计共有多少条reads（pair-end-reads这里算一条）参与了比对参考基因组

[sunchengquan 09:59:59 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |wc
  20000  391929 7049181

mate1 和 mate2 算一条read
$ cat tmp.sam |grep -v '^@' |cut -f1 |uniq |wc -l
10000

统计共有多少种比对的类型（即第二列数值有多少种）及其分布

[sunchengquan 10:10:04 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |cut -f2 |sort |uniq -c|sort -nrk1,1
   4650 83
   4650 163
   4516 99
   4516 147
    213 77
    213 141
    165 69
    165 137
    153 73
    153 133
    136 89
    136 165
    125 153
    125 101
     24 65
     24 129
     16 177
     16 113
      2 81
      2 161

筛选出比对失败的reads，看序列特征


[sunchengquan 10:19:05 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print}'|wc -l
1005
[sunchengquan 10:20:28 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print}'|cut -f10 |head -3
GCGTNGTACNNGNNGATGNTTACGACTAANACTTNNNNNTNCNGNNT
ANTGGNTNAGCGGNCGNAGNGNTCNGNAGCNANNAANCTGNTATGNTNNNTNNATCNNNCNNNNGTGNNAANNNCCAGNNNNGGCNNNNCATACANNNTNNNCNTNNNTACANNATANNCGNNATTNNACNGTNNCGNGNNCCGNGTNNNNN
NNTNNNNNAGNTNAANANANGTNTTGNGAANCTNNANACNNCTTTTNNNNTATGNTNNTNAGCCTAG
[sunchengquan 10:20:37 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$10}'|head -3
GCGTNGTACNNGNNGATGNTTACGACTAANACTTNNNNNTNCNGNNT
ANTGGNTNAGCGGNCGNAGNGNTCNGNAGCNANNAANCTGNTATGNTNNNTNNATCNNNCNNNNGTGNNAANNNCCAGNNNNGGCNNNNCATACANNNTNNNCNTNNNTACANNATANNCGNNATTNNACNGTNNCGNGNNCCGNGTNNNNN
NNTNNNNNAGNTNAANANANGTNTTGNGAANCTNNANACNNCTTTTNNNNTATGNTNNTNAGCCTAG

比对失败的reads区别成单端失败和双端失败情况，并且拿到序列ID

[sunchengquan 10:23:58 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort |uniq -c|grep -w 2 |head -2
      2 r1019
      2 r1105

[sunchengquan 10:24:20 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort |uniq -c|grep -wv 2 |head -2
      1 r1007
      1 r1010

筛选出比对质量值大于30的情况（看第五列）

[sunchengquan 10:24:42 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($5>30)print$1,$5}'|head -3
r1 42
r1 42
r2 42

筛选出比对成功，但是并不是完全匹配的序列

[sunchengquan 10:26:19 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($6!="*")print$6}'|grep '[IDNSHPX]'|head -4
72M2D119M
126M3D61M
70M1D71M
53M5D134M

筛选出insert size长度大于1250bp的pair-end reads

[sunchengquan 10:48:55 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cat tmp.sam |grep -v '^@' |awk '{if($9>1250||$9<-1250)print}'|less -S

统计参考基因组上面各条染色体的成功比对reads数量

[sunchengquan 22:25:45 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ cut -f3 tmp.sam |sort |uniq -c
    426 *
  19574 gi|9626243|ref|NC_001416.1|

测试数据只有一个染色体

筛选出原始fq序列里面有N的比对情况

[sunchengquan 11:35:18 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($10 ~/N/) print $0}' tmp.sam |less -S

筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况

[sunchengquan 11:42:42 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($10 ~/N/) print}' tmp.sam |awk '{if($6 !~ /[IDNSHPX]/) print}' |less -S
[sunchengquan 11:42:52 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($10 ~/N/) print}' tmp.sam |awk '{if($6 !~ /[IDNSHPX*]/) print}' |less -S

sam文件里面的头文件行数

[sunchengquan 12:37:28 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep '^@' tmp.sam 
@HD    VN:1.0    SO:unsorted
@SQ    SN:gi|9626243|ref|NC_001416.1|    LN:48502
@PG    ID:bowtie2    PN:bowtie2    VN:2.3.4.3    CL:"/home/sunchengquan/local/app/bowtie2-2.3.4.3-linux-x86_64/bowtie2-align-s --wrapper basic-0 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq"

sam文件里面每一行的tags个数一样吗？

[sunchengquan 13:43:16 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam |cut -f 12-1000|head -6
AS:i:-5    XN:i:0    XM:i:3    XO:i:0    XG:i:0    NM:i:3    MD:Z:59G13G21G26    YS:i:-4  YT:Z:CP
AS:i:-4    XN:i:0    XM:i:4    XO:i:0    XG:i:0    NM:i:4    MD:Z:65T3A19T107T0    YS:i:-5    YT:Z:CP
AS:i:-14    XN:i:0    XM:i:8    XO:i:0    XG:i:0    NM:i:8    MD:Z:0A0C0G0A108C23G9T81T46    YS:i:-5    YT:Z:CP
AS:i:-5    XN:i:0    XM:i:2    XO:i:0    XG:i:0    NM:i:2    MD:Z:1C182C3    YS:i:-14    YT:Z:CP
AS:i:-22    XN:i:0    XM:i:8    XO:i:0    XG:i:0    NM:i:8    MD:Z:80C4C16A52T23G30A8T76A41    YS:i:-3    YT:Z:CP
AS:i:-3    XN:i:0    XM:i:3    XO:i:0    XG:i:0    NM:i:3    MD:Z:23C1T39G1    YS:i:-22    YT:Z:CP

sam文件里每一行的tags个数分别是多少？

[sunchengquan 13:50:54 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam |cut -f 12-1000|awk '{print NF}' |sort |uniq -c
    457 1
    548 2
    579 8
  18416 9

[sunchengquan 13:59:09 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam |cut -f 12-1000|awk '{split($0,a, " ");print length(a)}' |sort|uniq -c
    457 1
    548 2
    579 8
  18416 9

sam文件里记录的参考基因组染色体长度分别是


[sunchengquan 14:02:12 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep  '^@SQ' tmp.sam |cut -f3|cut -d: -f2
48502

找到比对情况有insertion情况的


[sunchengquan 14:05:50 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6 ~/I/)print}' tmp.sam |less -S

至少有3个I
[sunchengquan 14:12:07 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6 ~/(I.*){3,}/)print}' tmp.sam |less -S

找到比对情况有deletion情况的

[sunchengquan 14:12:27 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6 ~/D/)print}' tmp.sam |less -S

取出位于参考基因组某区域的比对记录，比如5013到50130区域

[sunchengquan 14:13:04 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($4>5013 && $4<50130) print}' tmp.sam |less -S

把sam文件按照染色体及起始坐标排序

[sunchengquan 14:30:02 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ grep -v '^@' tmp.sam | awk '{print $4":" $0}' |sort -t: -nrk1,1|less -S

找到102M3D11M的比对情况，计算其reads片段长度

[sunchengquan 14:40:31 ~/local/app/bowtie2-2.3.4.3-linux-x86_64/example/reads]
$ awk '{if($6=="102M3D11M")print $1,"read length:"length($10)}' tmp.sam
r284 read length:113

sunchengquan

关注

3
点赞
踩
30

收藏

觉得还不错? 一键收藏
0
评论
生物信息数据格式：sam,bam格式

数据获取首先安装bowtie短序列比对软件wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zipunzip bowtie2-2.3.4.3-linux-x86_64.zipln -s ~/local/app/bowtie2-2.3.4.3...
复制链接

扫一扫

专栏目录