samtools使用说明

2018-12-08 18:18

名称: samtools Sequence Alignment/Map (SAM)格式的应用程序语法:

samtools view ‐bt ref_list.txt ‐o aln.bam aln.sam.gz samtools sort aln.bam aln.sorted samtools index aln.sorted.bam

samtools view aln.sorted.bam chr2:20,100,000‐20,200,000 samtools merge out.bam in1.bam in2.bam in3.bam samtools faidx ref.fasta

samtools pileup ‐f ref.fasta aln.sorted.bam samtools tview aln.sorted.bam ref.fasta

描述:

Samtools是一系列处理BAM格式序列的应用。它从SAM(Sequence Alignment/Map)格式输入或者输出为SAM格式，可以进行排序，合并和建立索引，并且允许快速地检索任意区域的读段（reads）。

Samtools工作在数据流中。它需要一个输入文件（-代表标准输入stdin，1>）和一个输出文件（-代表标准输出stdout）。因此，一些命令可以在unix管道结合。Samtools通常将警告信息和错误信息输出到标准错误输出（stderr，2>）。

Samtools也能打开远程FTP或者HTTP服务器上的BAM文件（BAM文件名起始为`ftp://'或者`http://'）。Samtools将会检查当前工作目录中的index文件，如果确实将会下载。除非要求它检索整个序列文件，否则不会这样做。

命令和选项：

1）samtools view *‐bhuHS+ *‐t in.refList+ *‐o output+ *‐f reqFlag+ *‐F skipFlag+ *‐q minMapQ+ *‐l library+ *‐r read‐Group+ | *region1 *...++，提取/打印所有的或者部分序列文件。如果没有指定区域，所有的序列都将会被打印。否则，仅仅打印覆盖特定区域的序列才会被打印出来。如果一个序列覆盖多个区域，那么它将会被多次打印出来。一个区域可以以如下例格式呈现：`chr2' (the whole chr2), `chr2:1000000' (region starting from 1,000,000bp) or `chr2:1,000,000‐2,000,000' (region between 1,000,000 and 2,000,000bp including the end points).（即染色体名：起始打印位点..终止打印位点）位置是以1开始记录第一个碱基的方式定位（1-based coordinate）。

-b 以BAM格式输出，可以用于samtools的后续分析

-u 以未压缩的BAM格式输出，可以节约时间，一般在管道执行时使用 -h 在结果中包含头header

-H 只输出头

-S 输入文件为SAM格式，如果确实@SQ头，则需要-t选项

-t FILE 这个文件是以TAB分隔的，每行必须包含参考序列名字（如染色体名），一行只包含单个参考序列名字，其它区域被忽略。这个文件定义排序是的参考序列顺序。如果你运行samtools faidx ，得到的索引结果文件可以作为这里的FILE文件（见6）。

-o FILE 输出文件[stdout]，如果在该管道中可以使用command1|tee FILE|command2。 -f INT 只输出在FLAG中含有整个INT的序列，INT可以使用十六进制的/^0x*0‐9A‐F++/ *默认0]

-F INT 跳过含有INT的序列

-q INT 跳过MAPQ（定位的质量）比INT小的序列[默认0] -l STR 只输出STR库中的序列

-r STR 只输出在STR读段组中的序列

2）samtools import ，从0.1.4起，它与samtools view ‐bt ‐o 意义相同。

3）samtools sort *‐n+ *‐m maxMem+ ，根据左起位点对序列排序，将产生.bam文件，当整个序列无法完全装载到内存（‐m选项控制）时，可能还会产生.%d.bam这个临时文件。

-n 设定排列方式为读段名称而非染色体位置 -m 设定内存使用量[500000000](默认500M)

4）samtools merge *‐h inh.sam+ *‐n+ [...]，将多个排序后的序列文件合并。参考头文件列出所有的输入BAM文件和inh.sam中@SQ头，如果存在，必须全部指向同一系列的参考序列。参考头文件列表（除非被-h选项替代）和in1.bam中的“@”头将会被拷贝到out.bam中，其他的输入文件中的头将会被忽略。

-h FILE 使用FILE中的行作为@头拷贝到out.bam中，代替从in1.bam中拷贝的头。（FILE

实际上是BAM格式的。但是其中的任何序列信息都会被忽略）。

-n 输入文件是以读段名称排序的，而非染色体定位（当sort中使用-n选项时，此处应该选-n）。

5）samtools index 对排序后的序列索引，用于快速的随机处理，此处将产生.bai文件。

6）samtools faidx [region1 [...]]，对FASTA格式的参考基因组序列建立索引或者提取已经索引的序列中的子序列（subsequence）。如果没有指定区域，它将会索引文件，并在磁盘上建立.fai，子序列将会被检索并且以FASTA格式打印在标准输出。输入文件可以压缩成RAZF格式。

7）samtools pileup *‐f in.ref.fasta+ *‐t in.ref_list+ *‐l in.site_list+ *‐iscgS2+ *‐T theta+ *‐N nHap+ *‐r pairDiffRate] |，以pileup格式输出序列，在pileup格式中，每一行代表一个基因组位点，包括染色体名称、位点、参考碱基、读段质量和序列定位质量。匹配、错配、插入和缺失、链、mapping质量以及读段的开始和结束位点都在读段碱基一列。在这一列中，点代表与参考序列的正向链匹配，逗号表示与参考序列的反向链匹配，ACTGN代表正向链上的错配，而actgn代表反向链上的错配，+[0‐9++*ACGTNacgtn++这一模式表示在参考序列之间存在插入序列，模式中，数字代表插入长度，其后为插入序列。‐*0‐9++*ACGTNacgtn++这一模式代表与参考序列相比存在缺失，后面以*代表被删除的碱基。在读段碱基列中，符号^标记一个读段片段的起始，它是读段中被N/S/H CIGAR操作符分隔开临近的子序列（subsequence）。^后面紧跟的ASCII码字符减去33给出mapping质量。$符号标记读段片段（read segment）的结尾。（注意-c选项不同的输出）

如果应用-c选项，一致碱基、Phred-scaled一致性质量、SNP质量（指与参考序列相同的

Phred-scaled一致性的可能性）和包含该位点的读段mapping质量均方根（root mean square RMS）将会被插入到参考碱基和读段序列之间（即会多出几列）。插入缺失将占有额外一列。每个插入缺失行包含染色体名称、位点、*、基因型、一致性质量、SNP质量、mapping质量均方根，#支持第一种等位序列的读段，#支持第二种等位序列的读段，#与前两种等位形式不同且包含插入缺失的读段。

-s 最后一行打印mapping质量。这一选项使得输出结果易于分析，但较耗费空间 -S 输入文件时SAM格式

-i 只输出含有插入缺失的pileup行

-f FILE FASTA格式的参考序列，如果缺少索引文件将产生FILE.fai -M INT设定mapping质量上限（cap）为INT

-t FILE以import 命令描述的格式列出参考名称和序列长度列表。如果出现这一个选项，samtools将假设输入文件为SAM格式；否则假设为BAM格式。

-l(L) FILE指定pileup输出的位点，第一列为染色体，第二列为1-based位点。其他的列将会被忽略，推荐与-s选项同时使用，因为默认格式我们可能无法知道mapping质量 -c 使用MAQ一致性模型检测一致序列，在使用-g和-c选项时，仅-T、-N、-l和-r选项可用 -g 产生GLFv3格式的基因型相似性，这一选项抑制-c、-i和-s选项 -T FLOAT MAQ一致性检测模式的theta参数[0.85]（错误依赖系数） -N INT样本中的单倍型数量[默认2，要求>=2]

-r FLOAT 每一对单倍体之间期望的片段差异值[0.001]

-I(i) INT 序列中插入缺失的Phred概率[40]

8）samtools tview [ref.fasta]，文本定位查看器，在查看器中，使用？取得帮助，按下g从一个区域开始位置检查定位（alignment），格式如chr10:10,000,000。 9）samtools fixmate ，为以名称排序的定位alignment填入配对坐标，ISIZE（inferred insert size猜测的插入序列大小）和配对相应的标签（flag）

10）samtools rmdup ，移除潜在的PCR重复：如果多个读段含有相同的外部位点（external coordinates），只保留那些高mapping质量的碱基对。这一命令只适用于FR开始的并且要求正确设置ISIZE。

11）samtools rmdupse ，除去单端读段（single-ended reads）。这一命令将把所有的读段当作单端读段，即使它们实际上是配对的。

12）samtools fillmd *‐e+ ，产生MD标签，如果已经存在并且产生的MD标签不同，这个命令将给予警告。

-e 如果读段碱基与参考序列相同，将其转换成=，indel检测目前还不支持=作为碱基标

志。

SAM格式：

SAM是以TAB分隔的，除去以@开始的header行，每一定位行（alignment）包含以下项目： |Col | Field | Description |

+‐‐‐‐+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+ | 1 | QNAME | Query (pair) NAME 读段对应模板名称| | 2 | FLAG | bitwise FLAG 按位的标签（意义见后面）|

| 3 | RNAME | Reference sequence NAME 参考序列名称|

| 4 | POS | 1‐based leftmost POSition/coordinate of clipped sequence 1-based左起位置/整齐的

序列定位|

| 5 | MAPQ | MAPping Quality (Phred‐scaled) 读段定位质量|

| 6 | CIAGR | extended CIGAR string 扩展的CIGAR字符串（详见SAM格式详解）|

| 7 | MRNM | Mate Reference sequence NaMe (`=' if same as RNAME) 配对的参考序列名称，=表示相同|

|10 | SEQ | query SEQuence on the same strand as the reference 与参考序列处于同一链的读段序列|

|11 | QUAL | query QUALity (ASCII‐33 gives the Phred base quality) 度短序列的碱基质量ASCII码-33为碱基质量|

|12 | OPT | variable OPTional fields in the format TAG:VTYPE:VALUE 变量选项，格式TAG:VTYPE:VALUE |

+‐‐‐‐+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+ 下表列出了FLAG列的含义：

+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+ | Flag | Description |

+‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+

|0x0001 (1)| the read is paired in sequencing |读段序列是成对的

|0x0100 (256)| the alignment is not primary |定位不是最优选

|0x0200 (512)| the read fails platform/vendor quality checks |读段质量未生成

|0x0400 (1024)| the read is either a PCR or an optical duplicate |读段是PCR或者光学重复 +‐‐‐‐‐‐‐+‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐+

限制：

1、在bam_import.c, bam_endian.h, bam.c和bam_aux.c中的非定位单词 2、 CIGAR操作符P不能正确处理

3、合并时，输入文件需要有相同数目的参考基因序列。设备要求可以降低，另外，合并时

不会自动重构头字典（header dictionary）。用户必须提供正确的头，Picard做合并表现更好。

4、 samtools的rmdup无法处理单端数据，不能移除染色体之间的重复。Picard表现更好。 AUTHOR

Heng Li from the Sanger Institute wrote the C version of samtools. Bob Handsaker from the Broad Institute implemented the BGZF library and Jue Ruan from Beijing Genomics Institute wrote the RAZF library. Various

people in the 1000Genomes Project contributed to the SAM format specification.