samtools 工具

参考文章:https://www.cnblogs.com/xiaofeiIDO/p/6805373.html
1、view
主要功能:sam和bam文件之间相互转换,针对bam文件进行相关操作。bam文件是sam文件的二进制格式,占据内存较小且运算速度快。
-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
7、flagstat

作用:reads的比对情况统计

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

用法:

samtools flagstat test.sort.bam
8、depth

作用:每个碱基位点的测序深度
-a:输出所有的碱基深度(包括0)

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

-f:bam文件名字

-l:设置read长度阈值

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

-q:碱基质量阈值

-Q:比对质量阈值

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

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、fastq

作用:bam文件转换为fastq
用法:

samtools fastq test.bam
12、fasta

作用:bam文件转换为fasta
用法:

samtools fasta test.bam
13、idxstats

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

Usage: samtools idxstats <in.bam>

用法:

samtools idxstats test.sort.bam

结果返回一个表格,4列。
第一列:序列名

第二列:序列长度

第三列:比对上的reads数

第四列:未必对数目

14、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的长度分布

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 202,607评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,047评论 2 379
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,496评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,405评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,400评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,479评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,883评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,535评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,743评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,544评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,612评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,309评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,881评论 3 306
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,891评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,136评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,783评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,316评论 2 342

推荐阅读更多精彩内容