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

相关文章

业务中14个需要进行A/B测试的时刻[信息图]

业务中14个需要进行A/B测试的时刻[信息图]

在本指南中，我们将全面了解有关 A/B测试 的所有内容。 我们将介绍不同类型的A/B测试，如何有效地规划和启动测试，如何评估测试是否成功，您应该关注哪些指标，多年来我们发现的常见错误等等。 什么是A/B测试？ A/B测试（有时称为“分割测试”）是一种实验类型，其中您创建两种或多种内容变体——如登录页面、电子邮件或广告——并将它们显示给不同的受众群体，以查看哪一种效果最好。 本质上，A/B测
阅读更多...
【北交大信息所AI-Max2】使用方法

【北交大信息所AI-Max2】使用方法

BJTU信息所集群AI_MAX2使用方法 使用的前提是预约到相应的算力卡，拥有登录权限的账号密码，一般为导师组共用一个。 有浏览器、ssh工具就可以。 1.新建集群Terminal 浏览器登陆10.126.62.75 （如果是1集群把75改成66） 交互式开发 执行器选Terminal 密码随便设一个（需记住） 工作空间：私有数据、全部文件 加速器选GeForce_RTX_2080_Ti
阅读更多...
easyui同时验证账户格式和ajax是否存在

easyui同时验证账户格式和ajax是否存在

accountName: {validator: function (value, param) {if (!/^[a-zA-Z][a-zA-Z0-9_]{3,15}$/i.test(value)) {$.fn.validatebox.defaults.rules.accountName.message = '账户名称不合法（字母开头，允许4-16字节，允许字母数字下划线）';return fal
阅读更多...
生信代码入门：从零开始掌握生物信息学编程技能

生信代码入门：从零开始掌握生物信息学编程技能

少走弯路,高效分析;了解生信云,访问 【生信圆桌x生信专用云服务器】 : www.tebteb.cc 介绍 生物信息学是一个高度跨学科的领域，结合了生物学、计算机科学和统计学。随着高通量测序技术的发展，海量的生物数据需要通过编程来进行处理和分析。因此，掌握生信编程技能，成为每一个生物信息学研究者的必备能力。 生信代码入门，旨在帮助初学者从零开始学习生物信息学中的编程基础。通过学习常用
阅读更多...
生信圆桌x生信分析平台：助力生物信息学研究的综合工具

生信圆桌x生信分析平台：助力生物信息学研究的综合工具

介绍 少走弯路,高效分析;了解生信云,访问 【生信圆桌x生信专用云服务器】 : www.tebteb.cc 生物信息学的迅速发展催生了众多生信分析平台，这些平台通过集成各种生物信息学工具和算法，极大地简化了数据处理和分析流程，使研究人员能够更高效地从海量生物数据中提取有价值的信息。这些平台通常具备友好的用户界面和强大的计算能力，支持不同类型的生物数据分析，如基因组、转录组、蛋白质组等。
阅读更多...
[数据集][目标检测]血细胞检测数据集VOC+YOLO格式2757张4类别

[数据集][目标检测]血细胞检测数据集VOC+YOLO格式2757张4类别

数据集格式：Pascal VOC格式+YOLO格式(不包含分割路径的txt文件，仅仅包含jpg图片以及对应的VOC格式xml文件和yolo格式txt文件) 图片数量(jpg文件个数)：2757 标注数量(xml文件个数)：2757 标注数量(txt文件个数)：2757 标注类别数：4 标注类别名称:["Platelets","RBC","WBC","sickle cell"] 每个类别标注的框数：
阅读更多...
Linux命令（11）：系统信息查看命令

Linux命令（11）：系统信息查看命令

系统 # uname -a # 查看内核/操作系统/CPU信息# head -n 1 /etc/issue # 查看操作系统版本# cat /proc/cpuinfo # 查看CPU信息# hostname # 查看计算机名# lspci -tv # 列出所有PCI设备# lsusb -tv
阅读更多...
一步一步将PlantUML类图导出为自定义格式的XMI文件

一步一步将PlantUML类图导出为自定义格式的XMI文件

一步一步将PlantUML类图导出为自定义格式的XMI文件 说明： 首次发表日期：2024-09-08PlantUML官网： https://plantuml.com/zh/PlantUML命令行文档： https://plantuml.com/zh/command-line#6a26f548831e6a8cPlantUML XMI文档: https://plantuml.com/zh/xmi
阅读更多...
【小迪安全笔记 V2022 】信息打点9~11

【小迪安全笔记 V2022 】信息打点9~11

第9天 信息打点-CDN绕过篇&漏洞回链8接口探针&全网扫指&反向件 知识点： 0、CDN知识-工作原理及阻碍 1、CDN配置-域名&区域&类型 2、CDN绕过-靠谱十余种技战法 3、CDN绑定-HOSTS绑定指向访问 CDN 是构建在数据网络上的一种分布式的内容分发网。 CDN的作用是采用流媒体服务器集群技术，克服单机系统输出带宽及并发能力不足的缺点，可极大提升系统支持的并发流数目，减少或避
阅读更多...
Weex入门教程之4，获取当前全局环境变量和配置信息(屏幕高度、宽度等)

Weex入门教程之4，获取当前全局环境变量和配置信息(屏幕高度、宽度等)

$getConfig() 获取当前全局环境变量和配置信息。 Returns: config (object): 配置对象；bundleUrl (string): bundle 的 url；debug (boolean): 是否是调试模式；env (object): 环境对象； weexVersion (string): Weex sdk 版本；appName (string): 应用名字；
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2