嫌picard去重复太慢?快来试试这两个高效快速的工具吧~
SAMBAMBA
http://lomereiter.github.io/sambamba/
一、安装 (支持mac OS/linux 64位)
git clone --recursive https://github.com/lomereiter/sambamba.git
cd sambamba
make
二、使用方法
1.排序
sambamba sort
OPTIONS <input.bam>
主要参数:
-o, --out 设置输出文件的名字(默认 .sorted.bam)
-n, --sort-by-name 按reads id排序(默认按照在参考基因组上的位置排序)
2.建立索引
sambamba index
[-p|--show-progress] [-n|--threads=NTHREADS] <input.bam> [<output.bai>]
示例:
$ sambamba index example.bam
#显示处理过程
$ sambamba index --show-progress example.bam /tmp/example.bam.bai
3.提取文件的信息
sambamba view
OPTIONS <input.bam | input.sam> [region1 [...]]
主要参数:
-S 输入文件为sam(默认为bam)
-F, --filter=FILTER 过滤提取bam
-f, --format=FORMAT 指定输出文件格式(默认是sam, 还支持bam, json, or msgpack )
-h, --with-header 保留header
示例:
#显示参考基因组序列基本信息
$ sambamba view --reference-info ex1_header.bam
[{"name":"chr1","length":1575},{"name":"chr2","length":1584}]
#计算3号染色体上质量值大于5且序列长大于80bp的reads个数
$ sambamba view -c -F "ref_id == 3 and mapping_quality >= 50 and sequence_length >= 80" ex1_header.bam
3124
4.合并多个bam文件
sambamba merge
OPTIONS <output.bam> <input1.sorted.bam> <input2.sorted.bam> [...] #automatically like in Picard merging tool
主要参数:
-t, --nthreads=NTHREADS 设置线程数
-H, --header 合并sam中的header
-l, --compression-level 按0 to 9设置文件压缩的程度
5.查看reads flag的比对结果
sambamba flagstat
OPTIONS <input.bam>
显示以下信息:
First line contains numbers of QC-passed and QC-failed reads. Then come pairs of numbers, the former for QC-passed reads, the latter for QC-failed ones:
- duplicates
- mapped reads (plus percentage relative to the numbers from the first line)
- reads with 'is_paired' flag set
- paired reads which are first mates
- paired reads which are second mates
- paired reads with 'proper_pair' flag set (plus percentage relative to the numbers of QC-passed/failed reads with 'is_paired' flag set)
- paired reads where both mates are mapped
- paired reads where read itself is unmapped but mate is mapped
- paired reads where mate is mapped to a different chromosome
- the same as previous but mapping quality is not less than
6.查重复序列
sambamba markdup
OPTIONS <input.bam> <output.bam>
主要参数:
-r, --remove-duplicates 去掉重复序列,否则仅mark重复序列
-t, --nthreads=NTHREADS
-l, --compression-level=N
--tmpdir=TMPDIR 设置中间文件生成目录(默认为/tmp)
此外,还可以提取sam文件的某一段,sambamba slice
OPTIONS <input.bam> region
SAMBLASTER
https://github.com/GregoryFaust/samblaster
https://doi.org/10.1093/bioinformatics/btu314
一、安装 (支持linux/mac OS Version 10.7以上)
git clone git://github.com/GregoryFaust/samblaster.git
cd samblaster
make
cp samblaster /usr/local/bin/.
二、使用方法
主要参数:
-i --input 输入sam文件(必须包含header且按reads id排序)
-o --output 输出sam文件
-d --discordantFile 输出discordant read pairs
-s --splitterFile 输出split reads
-u --unmappedFile 输出unmapped/clipped reads
其他参数:
-a --acceptDupMarks 不去重
-e --excludeDups 去掉discordant, splitter, and/or unmapped等重复(具体定义详见samblaster主页)
-r --removeDups 去掉重复(-e --excludeDups类似)
--addMateTags 添加MC and MQ tags
-M 与bwa mem -M 类似
示例:
#自动输出discordant read pairs和split read alignments:
bwa mem <idxbase> samp.r1.fq samp.r2.fq | samblaster -e -d samp.disc.sam -s samp.split.sam | samtools view -Sb - > samp.out.bam
#从bam文件中提取 split reads和discordants read pairs
samtools view -h samp.bam | samblaster -a -e -d samp.disc.sam -s samp.split.sam -o /dev/null
需要注意的是picard Markduplicates 和sambamba markdup的输入文件是bam格式,samblaster是sam格式