生物信息数据格式：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。

相关文章

HTML5表格语法格式详解

HTML5表格语法格式详解

《HTML5表格语法格式详解》在HTML语法中，表格主要通过table、tr和td3个标签构成,本文通过实例代码讲解HTML5表格语法格式，感兴趣的朋友一起看看吧... 目录一、表格1.表格语法格式2.表格属性 3.例子二、不规则表格1.跨行2.跨列3.例子一、表格在html语法中，表格主要通过< tab
阅读更多...
Linux下如何使用C++获取硬件信息

Linux下如何使用C++获取硬件信息

《Linux下如何使用C++获取硬件信息》这篇文章主要为大家详细介绍了如何使用C++实现获取CPU,主板,磁盘,BIOS信息等硬件信息,文中的示例代码讲解详细,感兴趣的小伙伴可以了解下... 目录方法获取CPU信息：读取"/proc/cpuinfo"文件获取磁盘信息：读取"/proc/diskstats"文
阅读更多...
Python将博客内容html导出为Markdown格式

Python将博客内容html导出为Markdown格式

《Python将博客内容html导出为Markdown格式》Python将博客内容html导出为Markdown格式,通过博客url地址抓取文章,分析并提取出文章标题和内容,将内容构建成html,再转... 目录一、为什么要搞？二、准备如何搞？三、说搞咱就搞！抓取文章提取内容构建html转存markdown
阅读更多...
如何自定义Nginx JSON日志格式配置

如何自定义Nginx JSON日志格式配置

《如何自定义NginxJSON日志格式配置》Nginx作为最流行的Web服务器之一,其灵活的日志配置能力允许我们根据需求定制日志格式,本文将详细介绍如何配置Nginx以JSON格式记录访问日志,这种... 目录前言为什么选择jsON格式日志？配置步骤详解1. 安装Nginx服务2. 自定义JSON日志格式各
阅读更多...
python dict转换成json格式的实现

python dict转换成json格式的实现

《pythondict转换成json格式的实现》本文主要介绍了pythondict转换成json格式的实现,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下... 一开始你变成字典格式data = [ { 'a' : 1, 'b' : 2, 'c编程' : 3,
阅读更多...
一文详解SQL Server如何跟踪自动统计信息更新

一文详解SQL Server如何跟踪自动统计信息更新

《一文详解SQLServer如何跟踪自动统计信息更新》SQLServer数据库中,我们都清楚统计信息对于优化器来说非常重要,所以本文就来和大家简单聊一聊SQLServer如何跟踪自动统计信息更新吧... SQL Server数据库中，我们都清楚统计信息对于优化器来说非常重要。一般情况下，我们会开启"自动更新
阅读更多...
Python如何获取域名的SSL证书信息和到期时间

Python如何获取域名的SSL证书信息和到期时间

《Python如何获取域名的SSL证书信息和到期时间》在当今互联网时代,SSL证书的重要性不言而喻,它不仅为用户提供了安全的连接,还能提高网站的搜索引擎排名,那我们怎么才能通过Python获取域名的S... 目录了解SSL证书的基本概念使用python库来抓取SSL证书信息安装必要的库编写获取SSL证书信息
阅读更多...
Win32下C++实现快速获取硬盘分区信息

Win32下C++实现快速获取硬盘分区信息

《Win32下C++实现快速获取硬盘分区信息》这篇文章主要为大家详细介绍了Win32下C++如何实现快速获取硬盘分区信息,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 实现代码CDiskDriveUtils.h#pragma once #include <wtypesbase
阅读更多...
Python中Windows和macOS文件路径格式不一致的解决方法

Python中Windows和macOS文件路径格式不一致的解决方法

《Python中Windows和macOS文件路径格式不一致的解决方法》在Python中,Windows和macOS的文件路径字符串格式不一致主要体现在路径分隔符上,这种差异可能导致跨平台代码在处理文... 目录方法 1：使用 os.path 模块方法 2：使用 pathlib 模块（推荐）方法 3：统一使
阅读更多...
Java中使用注解校验手机号格式的详细指南

Java中使用注解校验手机号格式的详细指南

《Java中使用注解校验手机号格式的详细指南》在现代的Web应用开发中,数据校验是一个非常重要的环节,本文将详细介绍如何在Java中使用注解对手机号格式进行校验,感兴趣的小伙伴可以了解下... 目录1. 引言2. 数据校验的重要性3. Java中的数据校验框架4. 使用注解校验手机号格式4.1 @NotBl
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2