Samtools 常用命令的总结

1. flags

1 0x1 这序列是PE双端测序
2 0x2 这序列和参考序列完全匹配,没有错配和缺失
4 0x4 这序列没有mapping到参考序列上
8 0x8 这序列的mate序列没有mapping到参考序列上
16 0x10 这序列比对到参考序列的负链上
32 0x20 这序列的mate序列比对到参考序列的负链上
64 0x40 这序列是read1
128 0x80 这序列是read2
256 0x100 这序列不是主要的比对,因为序列可能比对到参考序列的多个位置上
512 0x200 这序列没有通过QC
1024 0x400 这序列是PCR重复序列
2048 0x800 这序列是补充比对

2. view

samtools view [options] in.sam|in.bam|in.cram [region...]

-f 提取 ## -f 4 提取出没有mapping上的reads
-F 过滤 ## -F 4 过滤掉没有mapping上的reads,也就是说提取出mapping上的reads
-u 输出格式为未压缩的bam格式
-q 过滤掉MAPQ值低某个阈值 ## -q 1 过滤掉MAPQ值低于1的情况
-h 设定输出的SAM文件带有header
-b 输出格式设定为BAM
-S 输入格式为SAM

提取比对到参考序列的结果:

可以添加 -@ 24 开24线程,显著提高效率

samtools view -bF 4 tmp.bam > tmp_F.bam

提取双端序列都比对到参考序列(4+8)的结果:

samtools view -bF 12 tmp.bam > tmp_F.bam

提取比对到chr1的结果

samtools view -b tmp.bam chr1 > tmp_chr1.bam

注:With no options or regions specified, prints all alignments in the specified input alignment file (in SAM, BAM, or CRAM format) to standard output in SAM format (with no header),也就是说,没有设定输出格式的话,默认是输出SAM格式,并且是没有header的SAM

3. index

samtools index [-bc] [-m INT] aln.bam|aln.cram [out.index]

-b 创建一个bai索引,默认设定这个参数(如果在命令中没这个参数)

建索引(必须是已经使用默认排序后的):

samtools index tmp.bam

4. sort

samtools sort [-l level] [-m maxMem] [-o out.bam] [-O format] [-n] [-t tag] [-T tmpprefix] [-@ threads] [in.sam|in.bam|in.cram]

-m 设置内存使用大小,默认是500,000,000(现在支持K,M,G等缩写)
-n Sort by read names (i.e., the QNAME field) rather than by chromosomal coordinates(似乎一般也是使用默认排序,即Sort alignments by leftmost coordinates,因为index需要默认排序…)
-@ 设置线程数
-O 输出的格式(sam,bam,cram),默认是bam

使用默认排序:

sort -@ 5 tmp.bam >tmp.sorted.bam

5. merge

samtools merge [-nur1f] [-h inh.sam] [-R reg] [-b <list>] <out.bam> <in1.bam> [<in2.bam> <in3.bam> ... <inN.bam>]

-b 一个bam文件一行的一个bam list文件
-n The input alignments are sorted by read names rather than by chromosomal coordinates
-h Use the lines of FILE as `@’ headers to be copied to out.bam, replacing any header lines that would otherwise be copied from in1.bam. (FILE is actually in SAM format, though any alignment records it may contain are ignored.)
-c 当多个输入文件包含相同的@RG头ID时,只保留第一个到合并后输出的文件。当合并多个相同样本的不同文件时,非常有用
-p 与-c参数类似,对于要合并的每一个文件中的@PG ID只保留第一个文件中的@PG

merge前必须是已经sort的文件,如果只是单纯的merge:

samtools merge tmp1.bam tmp2.bam

6. mpileup

samtools mpileup [-EBugp] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]

从官方说明:Generate VCF, BCF or pileup for one or multiple BAM files可看出,可以用来和bcftools搭配Call SNPs

最常用的几个参数:

-f The faidx-indexed reference file in the FASTA format(有索引(faidx)文件的参考序列)
-l BED or position list file containing a list of regions or sites where pileup or BCF should be generated(bed格式的文件,如果需要只处理特定位点的bam文件的话)
-r Only generate pileup in region 搭配-l使用,比如可以指定染色体
-g Compute genotype likelihoods and output them in the binary call format (BCF).(输出bcf格式文件)
-u Generate uncompressed VCF/BCF output(如果后面接管道符的话,必须使用这个指定不进行压缩)

搭配bcftools使用:

samtools mpileup -ugf <ref.fa> <sample1.bam>| bcftools call -vmO z -o <study.vcf.gz>

7. tview

samtools tview [-p chr:pos] [-s STR] [-d display] <in.sorted.bam> [ref.fasta]

显示reads比对到基因组的情况,类似于基因组浏览器

8. faidx

samtools faidx <ref.fasta> [region1 [...]]

给参考序列建索引,或者从已建索引的参考序列中提取一定位置范围内的序列

9. depth

samtools depth [options] [in1.sam|in1.bam|in1.cram [in2.sam|in2.bam|in2.cram] [...]]

计算bam/sam文件每个位点的测序深度

10. flagstat

samtools flagstat in.sam|in.bam|in.cram

统计bam文件中reads的比对情况,如多少reads比对上等信息

samtools官网手册还介绍了其他好多的功能,具体可参见:
http://www.htslib.org/doc/samtools.html

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

作者:fqiang1024
链接:https://www.jianshu.com/p/dea4a6a59078
来源:简书
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

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