zoukankan      html  css  js  c++  java
  • [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

  • 相关阅读:
    开端
    springboot打包失败
    CONCAT_WS函数
    关于使用|作为分隔符
    JSONArray数组
    Math.ceil(double)向上取整
    $.unique(array)数组去重
    觉得没有问题,却始终没有按照预想的走的问题
    关于mouseover与mouseout以及mouseleave和mouseenter
    关于网页元素定义click事件,点击一次触发两次问题解决办法
  • 原文地址:https://www.cnblogs.com/xiaofeiIDO/p/6805373.html
Copyright © 2011-2022 走看看