sambamba与samtools的细节差异

生信菜鸟团

发布于 2022-05-24 16:31:44

3.7K10

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

3月份，我在生信菜鸟团的首次发文，假阳性突变的出现居然是因为duplicates mark的不够？，讲述了supplementary read不会被GATK MarkDuplicates标记为duplicates的问题。

之后，针对这个问题，我开始着手对手上的bam进行重处理，并写出通用流程供实验室使用。

流程

虽然上次我推荐了samtools rmdup和MarkDuplicatesSpark，但是考虑到大多数同学都更常使用GATK，而MarkDuplicatesSpark的速度实在是太慢，所以最终还是选择queryname排序后使用MarkDuplicates来处理。

将一个已有的bam重新mark duplicates的流程如下：

$ sambamba sort -t 24 -n tumor.bam -o query_sorted.bam
$ gatk MarkDuplicates -I query_sorted.bam -O gatk_markduplicates.bam -M temp.metrics
$ sambamba sort -t 24 gatk_markduplicates.bam -o gatk_markduplicates_sorted.bam
$ sambamba index gatk_markduplicates_sorted.bam

sambamba毕竟不如samtools常用，于是写了一个samtools版本。

$ samtools sort -@ 24 -n tumor.bam -o query_sorted.bam
$ gatk MarkDuplicates -I query_sorted.bam -O gatk_markduplicates.bam -M temp.metrics
$ samtools sort -@ 24 gatk_markduplicates.bam -o gatk_markduplicates_sorted.bam
$ samtools index gatk_markduplicates_sorted.bam

但是这个流程我之前没有测试过，测试一下。

居然出错了。

$ tail samtools_test.error
[Mon May 02 12:37:48 CST 2022] picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 25.35 minutes.
Runtime.totalMemory()=28222947328
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
java.lang.IllegalArgumentException: Alignments added out of order in SAMFileWriterImpl.addAlignment for file:///hpc/data/home/slst/zhangly/yangjk/duplicates/samtools/gatk_markduplicates.bam. Sort order is queryname. Offending records are at [A00583:225:H2FV5DSXY:1:1101:1000:9111] and [A00583:225:H2FV5DSXY:1:1101:1000:10207]
        at htsjdk.samtools.SAMFileWriterImpl.assertPresorted(SAMFileWriterImpl.java:197)
        at htsjdk.samtools.SAMFileWriterImpl.addAlignment(SAMFileWriterImpl.java:184)
        at htsjdk.samtools.AsyncSAMFileWriter.synchronouslyWrite(AsyncSAMFileWriter.java:36)
        at htsjdk.samtools.AsyncSAMFileWriter.synchronouslyWrite(AsyncSAMFileWriter.java:16)
        at htsjdk.samtools.util.AbstractAsyncWriter$WriterRunnable.run(AbstractAsyncWriter.java:123)
        at java.lang.Thread.run(Thread.java:745)

显示原因是queryname排序的问题，虽然排序是queryname，但是出现了冲突的record。

问题探索

既然用sambamba跑的时候成功了，说明sambamba的排序方式是GATK兼容的。

两个软件分别进行queryname排序生成结果。

$ sambamba sort -t 24 -n tumor.bam -o query_sorted_sambamba.bam
$ samtools sort -@ 24 -n tumor.bam -o query_sorted_samtools.bam

初步对比一下前面的queryname。

$ samtools view --no-header query_grouped_sambamba.bam | awk '{print $1}' | uniq | head -5
A00583:225:H2FV5DSXY:1:1101:10004:10363
A00583:225:H2FV5DSXY:1:1101:10004:11397
A00583:225:H2FV5DSXY:1:1101:10004:11584
A00583:225:H2FV5DSXY:1:1101:10004:14998
A00583:225:H2FV5DSXY:1:1101:10004:15217
$ samtools view --no-header query_grouped_samtools.bam | awk '{print $1}' | uniq | head -5
A00583:225:H2FV5DSXY:1:1101:1000:2722
A00583:225:H2FV5DSXY:1:1101:1000:2942
A00583:225:H2FV5DSXY:1:1101:1000:3881
A00583:225:H2FV5DSXY:1:1101:1000:4664
A00583:225:H2FV5DSXY:1:1101:1000:4946

前五项就有区别了，表明samtools和sambamba使用的是两种不同的排序方式。

sambamba中第一个是A00583:225:H2FV5DSXY:1:1101:10004:10363，samtools中第一个是A00583:225:H2FV5DSXY:1:1101:1000:2722。所以，初步分析，sambamba是认为4比:的ASCII码小，所以将同一位上带4的放在了前面。而samtools，显然是认为1000数值上比10004小，所以将1000放在了前面。这两种排序方式在使用sort命令的时候遇到过。

$ echo -e '200\n1000' | sort
1000
200
$ echo -e '200\n1000' | sort -n
200
1000

加了-n才按数值排序，默认是按ASCII码排序。

具体看一下错误信息中提到的queryname的前后内容。

$ samtools view --no-header query_grouped_sambamba.bam | awk '{print $1}' | uniq | grep -A 2 -B 2 -n -E "A00583:225:H2FV5DSXY:1:1101:1000:9111|A00583:225:H2FV5DSXY:1:1101:1000:10207"
20-A00583:225:H2FV5DSXY:1:1101:10004:5290
21-A00583:225:H2FV5DSXY:1:1101:10004:9706
22:A00583:225:H2FV5DSXY:1:1101:1000:10207
23-A00583:225:H2FV5DSXY:1:1101:1000:12054
24-A00583:225:H2FV5DSXY:1:1101:1000:12148
--
44-A00583:225:H2FV5DSXY:1:1101:1000:6605
45-A00583:225:H2FV5DSXY:1:1101:1000:6668
46:A00583:225:H2FV5DSXY:1:1101:1000:9111
47-A00583:225:H2FV5DSXY:1:1101:10013:10473
48-A00583:225:H2FV5DSXY:1:1101:10013:12759

在sambamba中，这两条reads隔得很远，而且排序顺序确实是按照每个位置的ASCII码排序的。

$ samtools view --no-header query_grouped_samtools.bam | awk '{print $1}' | uniq | grep -A 2 -B 2 -n -E "A00583:225:H2FV5DSXY:1:1101:1000:9111|A00583:225:H2FV5DSXY:1:1101:1000:10207"
8-A00583:225:H2FV5DSXY:1:1101:1000:6605
9-A00583:225:H2FV5DSXY:1:1101:1000:6668
10:A00583:225:H2FV5DSXY:1:1101:1000:9111
11:A00583:225:H2FV5DSXY:1:1101:1000:10207
12-A00583:225:H2FV5DSXY:1:1101:1000:12054
13-A00583:225:H2FV5DSXY:1:1101:1000:12148

samtools这里，两条reads连在一起，很明显是按照数值大小进行排序。这种方式虽然更直观，但是与GATK不兼容，所以GATK在看到之后就报错了。

试图解决

发现samtools的小问题之后，查阅了一下samtools-sort文档[1]。

“When the -n option is present, records are sorted by name. Names are compared so as to give a "natural" ordering — i.e. sections consisting of digits are compared numerically while all other sections are compared based on their binary representation. This means "a1" will come before "b1" and "a9" will come before "a10". Records with the same name will be ordered according to the values of the READ1 and READ2 flags (see flags).

原来官方已经写明了，samtools sort -n执行的排序方式是natural order，具体就是包含数字的部分会按照数值大小进行比较，导致a9在a10的前面。

我也在github搜了一下，看看有没有人跟我一样困惑为什么要这样做。发现真有人提这个问题，还不少。

比如smowton[2]就提出，samtools sort -n排序的结果不是lexical sort。lexical sort的意思是字典序，也就是按ASCII码，sambamba所采用的排序方式。官方的回答是看一下SAM/BAM sorting order clarification[3]。

经过我对这个issue的粗略阅读，发现：

官方早就知道这个问题。
官方没有认定哪个排序方式是真正正确的，所以不打算改。
官方选择的解决方案是在SAM格式的文档中阐述清楚这个问题（some definition of fixed）。

之后，我去看了一下提到的SAM格式文档[4]。

在Tag的SO部分。

“For queryname sort, no explicit requirement is made regarding the ordering other than that it be applied consistently throughout the entire file. Footnote 6: It is known that widely used software libraries have differing definitions of the queryname sort order, meaning care should be taken when operating on multiple files of varying provenance.

最后的建议是，不同软件有差异，在使用时要小心。

Section 1.3.1部分对常用的排序方式进行了介绍。

总结

今天遇到的问题其实并没有解决。总的来说，为了和GATK兼容，不应使用samtools进行queryname排序，要么使用sambamba，要么使用GATK SortSam（GATK自己肯定和自己兼容）。

最后，一个建议的对fastq处理的流程如下。sambamba中处理pipe需要指定/dev/stdin和/dev/stdout。

reference=~/project/0-reference/2-GATK/hg19/human_g1k_v37_decoy.fasta
name=SRR8619267
fq1=${name}_1.clean.fq.gz
fq2=${name}_2.clean.fq.gz
bwa mem -t 24 -R "@RG\tPU:"${name}"\tID:"${name}"\tSM:"${name}"\tLB:WXS\tPL:ILLUMINA" $reference $fq1 $fq2 | \
sambamba view -t 24 -f bam -S /dev/stdin -o /dev/stdout | sambamba sort -t 24 -n /dev/stdin -o ${name}.bam
gatk MarkDuplicates -I ${name}.bam -O ${name}_markdup.bam -M ${name}.metrics -ASO queryname
sambamba sort -t 24 ${name}_markdup.bam -o ${name}_sorted.bam
sambamba index ${name}_sorted.bam