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

相关文章

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
阅读更多...
Python批量调整Word文档中的字体、段落间距及格式

Python批量调整Word文档中的字体、段落间距及格式

《Python批量调整Word文档中的字体、段落间距及格式》这篇文章主要为大家详细介绍了如何使用Python的docx库来批量处理Word文档,包括设置首行缩进、字体、字号、行间距、段落对齐方式等,需... 目录关键代码一级标题设置  正文设置完整代码运行结果最近关于批处理格式的问题我查了很多资料，但是都没
阅读更多...
基于Python开发PDF转Doc格式小程序

基于Python开发PDF转Doc格式小程序

《基于Python开发PDF转Doc格式小程序》这篇文章主要为大家详细介绍了如何基于Python开发PDF转Doc格式小程序,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 用python实现PDF转Doc格式小程序以下是一个使用Python实现PDF转DOC格式的GUI程序，采用T
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2