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
来源:简书
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。
网友评论