[SAMtools] 常用指令总结

源自:http://sanwen.net/a/hirxmpo.html

samtools是一系列处理bam和sam格式文件的应用程序集合,具有众多的功能。

首先呢,bam和sam文件主要是bwa、bowtie、tophat等序列比对工具产生的,这些软件我们后面会谈到。

软件下载安装:

地址:https://sourceforge.net/projects/samtools/

解压下载后的压缩文件,然后你会看到README文件,里面有详细的安装操作说明。

安装成功后,运行samtools,你会看到:


目前最新版本是1.3.1

下面我们针对samtools的主要命令以及参数做个实例演示。

操作文件下载:

wget http://popgen.dk/software/download/angsd/bams.tar.gz

解压后,在bams文件夹下,你会看到10个bam文件:


名字太复杂,进行批量重命名

rename "s/.mapped.ILLUMINA.bwa.CEU.low_coverage.20111****14.bam//" *

结果如下:


1、view

主要功能:sam和bam文件之间相互转换,针对bam文件进行相关操作。bam文件是sam文件的二进制格式,占据内存较小且运算速度快。

查看view的主要参数:


重要参数释义:

-b:输出bam格式,用于后续分析

-C:输出CRAM文件

-1:快速压缩(需要-b)

-u:输出未压缩的bam文件,节约时间,占据较多磁盘空间(需要-b)

-h:默认输出sam文件不带表头,该参数设定后输出带表头信息

-H:仅仅输出表头信息

-c:仅打印匹配数

-o:输出文件(stdout标准输出)

-U:输出没有经过过滤选择的reads

-t:制表分隔符文件(需要提供额外的参考数据,比如参考基因组的索引)

-L:仅包括和bed文件有重叠的reads

-r:仅输出在STR读段组中的reads

-R:仅输出特定reads

-q:定位的质量大于INT[默认0]

-l:仅输出在STR 库中的reads

-F:获得比对上(mapped)的过滤设置[默认0]

-f:获得未必对上(unmapped)的过滤设置[默认0]

-T:使用fasta格式的参考序列

……

实例演示:

bam文件转换为sam文件

samtools view -h smallNA06985  > test.sam

sam文件转换为bam文件

samtools view -bS -1 test.sam > test.bam

提取比对到参考基因组上的数据

samtools view -bF 4 test.bam > test.F.bam

提取没有比对到参考基因组上的数据

samtools view -bf 4 test.bam > test.f.bam

双端reads都比对到参考基因组上的数据

samtools view -bF 12 test.bam > test.12.bam

单端reads1比对到参考基因组上的数据

samtools view -bF 4 -f 8  test .bam > test1.bam

单端reads2比对到参考基因组上的数据

samtools view -bF 8 -f 4 test.bam > test2.bam

 

2、sort

主要功能:对bam文件进行排序(不能对sam文件进行排序)


主要参数释义:

-l:设置文件压缩等级,0不压缩,9压缩最高

-m:每个线程运行内存大小(可使用K M G表示)

-n:按照read名称进行排序

-o:排序后的输出文件

-T:PREFIX临时文件前缀

-@:设置排序和压缩的线程数,默认单线程

用法:

samtools sort  -l 9 -m 90M -n -o test.sort.bam -T sorted -@ 2  test.bam

上述含义是:压缩最高级9、每一个线程内存90Mb、输出文件名test.sort.bam、临时文件前缀sorted、线程数2。

当然,最简单命令:

samtools sort test.bam -o test.sort.bam

 

3、index

主要功能:对bam文件建立索引,但在此之前必须进行排序(sort),生成后缀是.bai的文件。


参数释义:

-b:创建一个.bai格式的索引文件(默认)

-c:创建.csi格式的索引文件

-m:创建.csi文件,索引的最小间隔值

用法:

samtools index test.sort.bam

 

4、merge

功能:合并多个已经sort的bam文件

当有多个样本的bam文件时,可以使用samtools的merge命令将这些bam文件合并为一个排序的且保持所有输入记录并保持现有排序顺序的bam文件。


主要参数释义:

-n:输入根据read排序的文件

-r:RG标签添加到每个比对文件上,标签值来自文件名

-u:输出未压缩的bam文件

-f:覆盖同名文件

-1:压缩等级1

-l:压缩等级0-9

-R:合并输入文件的指定区域

-h:FILE 指定FILE内的’@’头复制到输出bam文件中并替换输出文件的文件头

-c:多个输入文件包含相同的@RG头ID时,只保留第一个到合并后输出的文件

-p:合并的每一个文件中的@PG ID只保留第一个文件中的@PG

-s:覆盖随机种子

-b:文件列表,一行一个

用法:

samtools merge merge.bam  smallNA06985.sort smallNA06994.sort

 

5、faidx

功能:对fasta格式的文件建立索引,后缀名.fai。根据索引文件和序列文件,可以快速提取任意区域的序列文件。

fasta序列格式要求:每条序列,除了最后一行外,其他行的长度必须相同!

为了方便,我们在NCBI上下载水稻NIP基因组的序列,进行演示:

地址:https://www.ncbi.nlm.nih.gov/genome/?term=rice

然后,进行解压缩,重命名为seuence.fa

用法:

samtools faidx sequence.fa

最后生成一个sequence.fa.fai索引文件,一共5列,每列之间tab分割。


第一列:序列的名称

第二列:序列长度

第三列:第一个碱基的偏移量,从0开始计数

第四列:除了最后一行外,序列中每行的碱基数

第五列:除了最后一行外,序列中每行的长度(包括换行符)

从中呢,我们可以有目的的提取序列:

提取水稻第一染色体:

samtools faidx sequence.fa Chr1 > Chr1.fa

提取水稻第一染色体100-200bp的序列:

samtools faidx sequence.fa Chr1:100-200 > Chr1_100_200.fa

 

6、tview

作用:直观显示reads比对到基因组的情况,和基因组浏览器有点类似。


-d:输出类型

-p:直接定位给定位置

-s:reads显示

当给出参考基因组的时候,会在第一排给出参考基因组的序列,否则第一排全用N表示。

首先利用sort进行排序后,在利用index建立索引后,用下面命令:

samtools tview test.sort.bam


具体操作:按下g键,如上界面,有一个Goto对话框,提示输入要到达基因组的某一个位点。在这里,由于没有提供参考基因组,第一排都是N。

我们可以使用H(左)J(上)K(下)L(右)移动显示界面。大写字母移动快,小写字母移动慢。使用空格键向左快速移动,CTRL+H向左移动1KB,CTRL+L向右移动1KB。

此外,可以利用颜色来标注比对质量、碱基质量、核苷酸等。

>=30的碱基质量或比对质量使用白色表示

20-39的用黄色表示

10-19的用绿色表示

0-9的用蓝色表示

使用字母r切换显示read name等。

具体的操作说明可以?键查看。

7、flagstat

作用:reads的比对情况统计

Usage: samtools flagstat [--input-fmt-option OPT=VAL] <in.bam>

用法:

samtools flagstat test.sort.bam



解释:

#总的reads

45456 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
321 + 0 duplicates

#总体reads比对率

45293 + 0 mapped (99.64% : N/A)

#PEreads
45327 + 0 paired in sequencing

#read1

22670 + 0 read1 

#read2

22657 + 0 read2 

#PE reads比对到同一条参考序列,并且两条reads之间的距离符合设置的阈值

44950 + 0 properly paired (99.17% : N/A) 

#PE中两条都比对到参考序列上的reads数

45001 + 0 with itself and mate mapped

#单独一条匹配到参考序列上的reads数,和上一个相加,则是总的匹配上的reads数。

163 + 0 singletons (0.36% : N/A)

#PE中两条分别比对到两条不同的参考序列的reads数

23 + 0 with mate mapped to a different chr

#PE中两条分别比对到两条不同的参考序列的reads数(比对质量>=5)

11 + 0 with mate mapped to a different chr (mapQ>=5)

8、depth

作用:每个碱基位点的测序深度


-a:输出所有的碱基深度(包括0)

-b/-r:控制深度的范围(后面跟染色体)

-f:bam文件名字

-l:设置read长度阈值

-d/-m:最大覆盖深度

-q:碱基质量阈值

-Q:比对质量阈值

samtools depth -a -r 3 test.sort.bam 


结果是tab分割的三列

第一列:染色体名称

第二列:碱基位置

第三列:测序深度

9、mpileup

作用:对参考基因组每个位点做碱基堆积,用于call SNP和INDEL。主要是生成BCF、VCF文件或者pileup一个或多个bam文件。比对记录以在@RG中的样本名作为区分标识符。如果样本标识符缺失,那么每一个输入文件则视为一个样本。


主要参数释义:

-A:在检测变异中,不忽略异常的reads对

-C:用于调节比对质量的系数,如果reads中含有过多的错配,不能设置为零

-D:输出每个样本的reads深度

-l:BED文件或者包含区域位点的位置列表文件

注意:位置文件包含两列,染色体和位置,从1开始计数。BED文件至少包含3列,染色体、起始和终止位置,开始端从0开始计数。

-r:在指定区域产生pileup,需已建立索引的bam文件,通常和-l参数一起使用

-o/g/v:输出文件类型(标准格式文件或者VCF、BCF文件)

-t:设置FORMAT和INFO的列表内容,以逗号分割

-u:生成未压缩的VCF和BCF文件

-I:跳过INDEL检测

-m:候选INDEL的最小间隔的reads

-f:输入有索引文件的fasta参考序列

-F :含有间隔reads的最小片段

……

用法:

生成一个简单的vcf文件

samtools mpileup -vu test.sort.bam

如果有参考基因组的话

samtools mpileup -vuf genome.fasta  test.sort.bam

 

 

10、dict

作用:建立参考基因组字典

用法:

samtools dict  test.sort.bam sequences.fa 

11、cat

作用:连接多个bam文件(不做排序)

Usage: samtools cat [-h header.sam] [-o out.bam] <in1.bam> [...]

用法:

samtools cat  -o out.bam smallNA06985 smallNA06994

12、split

作用:根据read group 分割bam文件


用法:

samtools split test.sort.bam 
13、quickcheck

作用:检查SAM/BAM/CRAM文件的完整性


用法:

samtools quickcheck -v *.bam > bad_bams

14、fastq

作用:bam文件转换为fastq


用法:

samtools fastq test.bam

15、fasta

作用:bam文件转换为fasta


用法:

samtools fasta test.bam

16、idxstats

作用:检索和打印与输入文件相对应的index file里的统计信息

Usage: samtools idxstats <in.bam>

用法:

samtools idxstats test.sort.bam

结果返回一个表格,4列。


第一列:序列名

第二列:序列长度

第三列:比对上的reads数

第四列:未必对数目

 

17、stats

作用:对bam文件做详细统计,其统计结果可用mics/plot-bamstats作图

用法:

samtools stats test.bam

输出的信息比较多,部分如下:
Summary Numbers,raw total sequences,filtered sequences, reads mapped, reads mapped and paired,reads properly paired等信息
Fragment Qualitites:根据cycle统计每个位点上的碱基质量分布
Coverage distribution:深度为1,2,3,,,的碱基数目
ACGT content per cycle:ACGT在每个cycle中的比例
Insert sizes:插入长度的统计
Read lengths:read的长度分布

18、reheader

作用:替换bam文件的头


用法:

samtools view -H test.bam > in.header.bam

samtools reheader -p -i in.header.bam test.bam

19、rmdup

作用:将由PCR duplicates 获得的reads去掉,并保留高比对质量的reads


用法:

samtools rmdup -sS test.bam  output.bam

 

20、phase

作用:call杂合SNP,确定相位


用法:

samtools phase test.bam

 

21、calmad

作用:计算MD tag(a optional field,记录了mismatch信息)


用法:

samtools calmd test.bam sequences.fa

转载于:https://www.cnblogs.com/xiaofeiIDO/p/6805373.html

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://xiahunao.cn/news/1381026.html

如若内容造成侵权/违法违规/事实不符,请联系瞎胡闹网进行投诉反馈,一经查实,立即删除!

相关文章

Glide 的超时控制相关处理

作者&#xff1a;newki 前言 Glide 相信大家都不陌生&#xff0c;各种源码分析&#xff0c;使用介绍大家应该都是烂熟于心。但是设置 Glide 的超时问题大家遇到过没有。 我遇到了&#xff0c;并且掉坑里了&#xff0c;情况是这样的。 调用接口从网络拉取用户头像&#xff0c…

高通量DNA测序数据的生物信息学方法

来源&#xff1a;大数据期刊 时间&#xff1a;2016-05-13 14:41:09 作者&#xff1a;詹晓娟 姚登举 朱怀球 詹晓娟1&#xff0c;姚登举2&#xff0c;朱怀球3 1. 黑龙江工程学院计算机科学与技术学院&#xff0c;黑龙江 哈尔滨 150050&#xff1b; 2. 哈尔滨理工大学软件学院…

C++ pair详解

pair pair在cplusplus 与CPrimer中的介绍 首先可以看到pair是一个class template —类模板 pair也是一种模板类型。 对 pair 的介绍是&#xff1a; 这个类将一对值耦合在一起&#xff0c;这些值可能是不同类型的(T1和T2)。单个值可以通过其公共成员first和second访问。 pair是…

聊聊OpenStack运维架构

前言 想一想&#xff0c;从事OpenStack杂七杂八的事儿&#xff0c;至今正好三年半了。做过QA测试&#xff08;手动的、自动的&#xff09;、CI&#xff08;gerrit、jenkins、gitlab、harbor&#xff09;、云产品封装&#xff08;从系统pxe到openstack代码&#xff09;、自动化…

Lua pairs与ipairs效率分析

介于大家目前有些人比较关心 lua table中pairs 和 ipairs的效率问题, 特此研究了一下... 如有不正 还需指出.. 首先来看下 lua中table的结构定义: table中分为2个存储空间, 一个是线性数组空间(TValue *array), 和一个hash空间(Node *node) 当我们使用 pairs 和 ipairs 会产生…

【机器学习|数学基础】Mathematics for Machine Learning系列之矩阵理论(22):方阵函数在微分方程组中的应用

目录 前言往期文章5.6 方阵函数在微分方程组中的应用5.6.1 解一阶线性常系数齐次微分方程组5.6.2 解一阶线性常系数非齐次微分方程组 结语 前言 Hello&#xff01;小伙伴&#xff01; 非常感谢您阅读海轰的文章&#xff0c;倘若文中有错误的地方&#xff0c;欢迎您指出&#xf…

C++中pair使用详细说明

一、pair 的介绍 pair 是一个很实用的 "小玩意"&#xff0c;当想要将两个元素绑在一起作为一个合成元素、又不想要因此定义结构体时&#xff0c;使用 pair 可以很方便地作为一个代替品。 也就是说&#xff0c;pair 实际上可以看作一个内部有两个元素的结构体&#xf…

Solidity实现默克尔树 Merkle Tree

​​Merkle Tree​​​&#xff0c;也叫默克尔树或哈希树&#xff0c;是区块链的底层加密技术&#xff0c;被BTC和Ethereum区块链广泛采用。​​Merkle Tree​​​是一种自下而上构建的加密树&#xff0c;每个叶子是对应数据的哈希&#xff0c;而每个非叶子为它的​​2​​个子…

论文浅尝 | PairRE: 通过成对的关系向量实现知识图谱嵌入

笔记整理&#xff1a;黎洲波&#xff0c;浙江大学硕士&#xff0c;研究方向为自然语言处理、知识图谱。 研究背景 知识图谱因其在问答、语义解析和命名实体消歧等任务取得了良好的效果而受到广泛关注&#xff0c;而大部分知识图谱都存在不全和缺失实体链接的问题&#xff0c;所…

致敬乔布斯的经典,锤子坚果Pro成2017年最受欢迎手机看罗永浩怎么说

锤子坚果Pro发布已经近2个月&#xff0c;但热度依旧不减。在刚刚过去的京东618活动中&#xff0c;坚果Pro在1500到2000元档位产品中一举斩获单品销量冠军。坚果Pro凭借出色的销售战绩坐实2017手机圈“黑马”之名&#xff0c;而其销量节节攀升&#xff0c;这其中必有一番原因。日…

坚果Pro 2安抚了不少人锤粉, 但用户更期待锤子T3

今年秋季&#xff0c;锤子科技创始人罗永浩于2017年11月7日在成都发布坚果系列2代手机坚果Pro 2&#xff0c;指纹和logo的融合增强了手机的一体型。软件上的再度优化&#xff0c;帮助盲人更注重人文关怀&#xff0c;再从罗永浩自带“流量”&#xff0c;坚果Pro2自然而然受到大家…

pro坚果android耗流量,深度使用坚果Pro3一个月,憋了一肚子话,不吐不快​

原标题&#xff1a;深度使用坚果Pro3一个月&#xff0c;憋了一肚子话&#xff0c;不吐不快​ 罗永浩创办的锤子手机曾经在国内手机市场&#xff0c;也是一枚耀眼的新兴&#xff0c;罗永浩对于手机工业设计的高标准严要求让锤子手机成为了国内少有的在设计上能和苹果三星比肩的手…

厉害了!原来这些文艺明星都喜欢锤子坚果Pro

最近在手机圈出现了一匹黑马&#xff0c;那就是锤子坚果Pro。在京东618期间取得了十分骄人的战果&#xff0c;荣获6月1日至18日1500-1999元价位档单品销量第一&#xff0c;成为该价格区间最受欢迎的手机&#xff0c;同时在2017年4月1日后首发的新品销量排名中位列第三&#xff…

锤子t1android驱动,锤子T1痛失安卓5.1!都是因为这?

现如今很多高端手机都开始升级安卓6.0了&#xff0c;但是情怀锤子却突然给了老用户一“锤子”&#xff0c;宣布第一代T1将不会升级到安卓5.1&#xff0c;因为“优化效果不明显”。 这顿时引发了一片争议。有的T1用户表示了理解&#xff0c;称手机够用就好&#xff0c;不在乎系统…

内蒙古大学计算机考研资料汇总

内蒙古大学研招网 内蒙古大学计算机学院 内蒙古大学计算机学院成立于1997年&#xff0c;其前身是1978年设置的计算机专业和1988年成立的计算机科学系。内蒙古大学软件学院成立于2005年&#xff0c;与计算机学院为一个实体&#xff0c;两个牌子。目前学院由计算机科学系、…

如何获取bainu文档并用斡仑office进行编码转换-永中office蒙文版

声明&#xff1a; 1.bainu软件是由内蒙古卓嘎信息技术有限公司研发的。 2.斡仑office是由内蒙古斡仑科技有限公司与永中软件股份有限公司联合开发的蒙汉多文种跨平台办公套装。 第一&#xff0c;首先我们打开bainu软件&#xff0c;如图&#xff1a; 第二&#xff0c;下列图中…

为什么 Mixin 被认为是有害的

为什么 Mixin 被认为是有害的 Mixin 是在 Vue 2 中引入的&#xff0c;作为组件之间共享代码的解决方案&#xff0c;这种方式成为许多代码库不可或缺的一部分。然而&#xff0c;随着时间的推移&#xff0c;它们的使用开始出现问题。尽管 mixins 最初很有吸引力&#xff0c;但现…

【从零开始学习JAVA | 第四十五篇】反射

目录 前言&#xff1a; ​反射&#xff1a; 使用反射的步骤&#xff1a; 1.获取阶段&#xff1a; 2.使用阶段&#xff1a; 反射的应用场景&#xff1a; 使用反射的优缺点&#xff1a; 总结&#xff1a; 前言&#xff1a; Java中的反射是一项强大而灵活的功能&#xff0…

1688采源宝的商家靠谱吗 怎么入驻成为阿里采源宝商家

说到采源宝&#xff0c;我想很多微商都是不陌生的&#xff0c;采源宝的主要作用就是方便微商查看并转发供应商所发布的商品&#xff0c;并在有客户下单时&#xff0c;还可以轻松向供货商去下单。但也有很多朋友对采源宝的商家靠谱吗这个问题存在很大疑惑&#xff0c;下面我们就…

开发nft数字藏品平台合法合规吗?

开发nft数字产品平台合法合规吗&#xff1f;这是很多人都在问的问题&#xff0c;那么今天就来给大家说说。 开发nft数字产品平台可以肯定的是合法合规。 其实很多人对国内国外的玩法不太了解&#xff0c;有很多人对数字产品法律和政策不太清楚。首先&#xff0c;我们要了解国内…