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

2024-02-15 12:58
文章标签 格式 数据格式 信息 sam 生物 bam

本文主要是介绍生物信息数据格式：sam,bam格式，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

文章目录

  • 数据获取
  • 格式说明
  • 实例演练
    • 统计共有多少条reads（pair-end-reads这里算一条）参与了比对参考基因组
    • 统计共有多少种比对的类型（即第二列数值有多少种）及其分布
    • 筛选出比对失败的reads，看序列特征
    • 比对失败的reads区别成单端失败和双端失败情况，并且拿到序列ID
    • 筛选出比对质量值大于30的情况（看第五列）
    • 筛选出比对成功，但是并不是完全匹配的序列
    • 筛选出insert size长度大于1250bp的pair-end reads
    • 统计参考基因组上面各条染色体的成功比对reads数量
    • 筛选出原始fq序列里面有N的比对情况
    • 筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况
    • sam文件里面的头文件行数
    • sam文件里面每一行的tags个数一样吗？
    • sam文件里每一行的tags个数分别是多少？
    • sam文件里记录的参考基因组染色体长度分别是
    • 找到比对情况有insertion情况的
    • 找到比对情况有deletion情况的
    • 取出位于参考基因组某区域的比对记录，比如5013到50130区域
    • 把sam文件按照染色体及起始坐标排序
    • 找到102M3D11M的比对情况，计算其reads片段长度

数据获取

首先安装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.zipunzip 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.sambowtie2 -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 '^@' |wc20000  391929 7049181mate1 和 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,14650 834650 1634516 994516 147213 77213 141165 69165 137153 73153 133136 89136 165125 153125 10124 6524 12916 17716 1132 812 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 -22 r10192 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 -21 r10071 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 -c426 *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 -c457 1548 2579 818416 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 -c457 1548 2579 818416 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

这篇关于生物信息数据格式：sam,bam格式的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/711453。

相关文章

C++ 函数 strftime 和时间格式示例详解

C++ 函数 strftime 和时间格式示例详解

《C++函数strftime和时间格式示例详解》strftime是C/C++标准库中用于格式化日期和时间的函数,定义在ctime头文件中,它将tm结构体中的时间信息转换为指定格式的字符串,是处理... 目录C++ 函数 strftipythonme 详解一、函数原型二、功能描述三、格式字符串说明四、返回值五
阅读更多...
C#实现将Office文档(Word/Excel/PDF/PPT)转为Markdown格式

C#实现将Office文档(Word/Excel/PDF/PPT)转为Markdown格式

《C#实现将Office文档(Word/Excel/PDF/PPT)转为Markdown格式》Markdown凭借简洁的语法、优良的可读性,以及对版本控制系统的高度兼容性,逐渐成为最受欢迎的文档格式... 目录为什么要将文档转换为 Markdown 格式使用工具将 Word 文档转换为 Markdown（.
阅读更多...
Java中JSON格式反序列化为Map且保证存取顺序一致的问题

Java中JSON格式反序列化为Map且保证存取顺序一致的问题

《Java中JSON格式反序列化为Map且保证存取顺序一致的问题》：本文主要介绍Java中JSON格式反序列化为Map且保证存取顺序一致的问题,具有很好的参考价值,希望对大家有所帮助,如有错误或未... 目录背景问题解决方法总结背景做项目涉及两个微服务之间传数据时，需要提供方将Map类型的数据序列化为co
阅读更多...
Linux查看系统盘和SSD盘的容量、型号及挂载信息的方法

Linux查看系统盘和SSD盘的容量、型号及挂载信息的方法

《Linux查看系统盘和SSD盘的容量、型号及挂载信息的方法》在Linux系统中,管理磁盘设备和分区是日常运维工作的重要部分,而lsblk命令是一个强大的工具,它用于列出系统中的块设备（blockde... 目录1. 查看所有磁盘的物理信息方法 1：使用 lsblk（推荐）方法 2：使用 fdisk -l（
阅读更多...
SpringBoot如何对密码等敏感信息进行脱敏处理

SpringBoot如何对密码等敏感信息进行脱敏处理

《SpringBoot如何对密码等敏感信息进行脱敏处理》这篇文章主要为大家详细介绍了SpringBoot对密码等敏感信息进行脱敏处理的几个常用方法,文中的示例代码讲解详细,感兴趣的小伙伴可以了解下... 目录​1. 配置文件敏感信息脱敏​​2. 日志脱敏​​3. API响应脱敏​​4. 其他注意事项​​总结
阅读更多...
Ubuntu上手动安装Go环境并解决“可执行文件格式错误”问题

Ubuntu上手动安装Go环境并解决“可执行文件格式错误”问题

《Ubuntu上手动安装Go环境并解决“可执行文件格式错误”问题》：本文主要介绍Ubuntu上手动安装Go环境并解决“可执行文件格式错误”问题,具有很好的参考价值,希望对大家有所帮助,如有错误或未... 目录一、前言二、系统架构检测三、卸载旧版 Go四、下载并安装正确版本五、配置环境变量六、验证安装七、常见
阅读更多...
springboot实现配置文件关键信息加解密

springboot实现配置文件关键信息加解密

《springboot实现配置文件关键信息加解密》在项目配置文件中常常会配置如数据库连接信息,redis连接信息等,连接密码明文配置在配置文件中会很不安全,所以本文就来聊聊如何使用springboot... 目录前言方案实践1、第一种方案2、第二种方案前言在项目配置文件中常常会配置如数据库连接信息、Red
阅读更多...
使用Python开发Markdown兼容公式格式转换工具

使用Python开发Markdown兼容公式格式转换工具

《使用Python开发Markdown兼容公式格式转换工具》在技术写作中我们经常遇到公式格式问题,例如MathML无法显示,LaTeX格式错乱等,所以本文我们将使用Python开发Markdown兼容... 目录一、工具背景二、环境配置（Windows 10/11）1. 创建conda环境2. 获取XSLT
阅读更多...
Go语言开发实现查询IP信息的MCP服务器

Go语言开发实现查询IP信息的MCP服务器

《Go语言开发实现查询IP信息的MCP服务器》随着MCP的快速普及和广泛应用,MCP服务器也层出不穷,本文将详细介绍如何在Go语言中使用go-mcp库来开发一个查询IP信息的MCP... 目录前言mcp-ip-geo 服务器目录结构说明查询 IP 信息功能实现工具实现工具管理查询单个 IP 信息工具的实现服
阅读更多...
使用Python从PPT文档中提取图片和图片信息(如坐标、宽度和高度等)

使用Python从PPT文档中提取图片和图片信息(如坐标、宽度和高度等)

《使用Python从PPT文档中提取图片和图片信息(如坐标、宽度和高度等)》PPT是一种高效的信息展示工具,广泛应用于教育、商务和设计等多个领域,PPT文档中常常包含丰富的图片内容,这些图片不仅提升了... 目录一、引言二、环境与工具三、python 提取PPT背景图片3.1 提取幻灯片背景图片3.2 提取
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2