bamRemoveDuplicates fails to call duplicates

[ME]

Hi there,

I have generated the following test file (using wigsim to simulate reads and STAR to align them to the mouse genome) which contains some duplicate reads (i.e. _2d3 and _8b5;  _52b, _45c and _432; _932 and _8bc; file is attached as test.bam):

@HD    VN:1.4    SO:coordinate
@SQ    SN:6    LN:149736546
@PG    ID:STAR    PN:STAR    VN:STAR_2.5.2b    CL:STAR   --runThreadN 3   --genomeDir /home/me/indices/STAR/reference.43   --readFilesIn simulated_reads1.fq      --outSAMtype BAM   SortedByCoordinate      --outFilterType BySJout
@CO    user command line: STAR --genomeDir /home/me/indices/STAR/reference.43 --readFilesIn simulated_reads1.fq --outFilterType BySJout --runThreadN 3 --outSAMtype BAM SortedByCoordinate
chr6:59352460-59426290(-)_2612_3099_2:0:0_3:0:0_a61    0    6    59352462    255    100M    *    0    0    AATATTTTCATTATATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACCTTAAGGTTCTTAAAAGCTACTTAAACATATGTGGTTC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:92    nM:i:3
chr6:59352460-59426290(-)_2653_3098_2:0:0_3:0:0_a74    0    6    59352463    255    100M    *    0    0    ATATTTTCATTTTATCTTATTTGTCCCTGGGATAACCTACATTTGTTTTTCAAGATTTCACATTAAGGTTCATAAAAGCTACTTAAACATTTGTGGTTCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:92    nM:i:3
chr6:59352460-59426290(-)_2552_3097_3:0:0_2:0:0_875    0    6    59352464    255    100M    *    0    0    TATTTTCATTTTATCTTATTTGTCCCTGGGATAACATAGATTTGTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTTTGGTTCTA    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:94    nM:i:2
chr6:59352460-59426290(-)_2601_3093_0:0:0_2:0:0_31    0    6    59352468    255    100M    *    0    0    TTCATTTTATCTTATTTGTGCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATTCAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCAT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:94    nM:i:2
chr6:59352460-59426290(-)_2668_3091_0:0:0_4:0:0_8fe    0    6    59352470    255    100M    *    0    0    CATTTTATCTTATTAGTCCCTGGGATAACATCCATTTGTTTTTCAAGTTTTCCCATTAAGGTTCTTCAAAGCTACTTAAACATTTGTGGTTCTAGCATCC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:90    nM:i:4
chr6:59352460-59426290(-)_2578_3090_1:0:0_2:0:0_2d3    0    6    59352471    255    100M    *    0    0    ATTTTATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAAGATTTGTGTTTCTAGCATCCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:94    nM:i:2
chr6:59352460-59426290(-)_2624_3090_1:0:0_2:0:0_8b5    0    6    59352471    255    100M    *    0    0    ATTTTATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTGACATTAAGGTTCTTAAAAGGTACTTAAACATTTGTGGTTCTAGCATCCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:94    nM:i:2
chr6:59352460-59426290(-)_2644_3089_3:0:0_3:0:0_52b    0    6    59352472    255    100M    *    0    0    TTTTATCTTATTTGTCCCTGGGCTAACATACATTAGATTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:92    nM:i:3
chr6:59352460-59426290(-)_2491_3089_0:0:0_3:0:0_45c    0    6    59352472    255    99M1S    *    0    0    TTTTATCTTATTTGTCCCTGGGATCACATACATTTGTTTTTCAAGTTTTCACATTAAGGTTGTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTG    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:93    nM:i:2
chr6:59352460-59426290(-)_2578_3089_1:0:0_3:0:0_432    0    6    59352472    255    100M    *    0    0    TTTAATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACAATAAGGTTCTTACAAGCTACTTAAACATTTGTGGTTCTAGCATCCTC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:92    nM:i:3
chr6:59352460-59426290(-)_2577_3088_0:0:0_1:0:0_648    0    6    59352473    255    100M    *    0    0    TTTATCTTATTTGTCCCTGGGATAACATACATTTTTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:96    nM:i:1
chr6:59352460-59426290(-)_2628_3086_4:0:0_2:0:0_982    0    6    59352475    255    100M    *    0    0    TATGTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATAAAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTCTTG    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:94    nM:i:2
chr6:59352460-59426290(-)_2716_3086_1:0:0_3:0:0_8bc    0    6    59352475    255    100M    *    0    0    TATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTGAGGTTGTAGCATCCACTTG    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1    HI:i:1    AS:i:92    nM:i:3

When I run

$ STAR --runMode inputAlignmentsFromBAM --inputBAMfile test.bam --bamRemoveDuplicatesType UniqueIdentical
$ samtools view Processed.out.bam

I get

chr6:59352460-59426290(-)_2612_3099_2:0:0_3:0:0_a61     0       6       59352462        255     100M    *       0       0       AATATTTTCATTATATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACCTTAAGGTTCTTAAAAGCTACTTAAACATATGTGGTTC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:92 nM:i:3
chr6:59352460-59426290(-)_2653_3098_2:0:0_3:0:0_a74     0       6       59352463        255     100M    *       0       0       ATATTTTCATTTTATCTTATTTGTCCCTGGGATAACCTACATTTGTTTTTCAAGATTTCACATTAAGGTTCATAAAAGCTACTTAAACATTTGTGGTTCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:92 nM:i:3
chr6:59352460-59426290(-)_2552_3097_3:0:0_2:0:0_875     0       6       59352464        255     100M    *       0       0       TATTTTCATTTTATCTTATTTGTCCCTGGGATAACATAGATTTGTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTTTGGTTCTA    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:94 nM:i:2
chr6:59352460-59426290(-)_2601_3093_0:0:0_2:0:0_31      0       6       59352468        255     100M    *       0       0       TTCATTTTATCTTATTTGTGCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATTCAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCAT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:94 nM:i:2
chr6:59352460-59426290(-)_2668_3091_0:0:0_4:0:0_8fe     0       6       59352470        255     100M    *       0       0       CATTTTATCTTATTAGTCCCTGGGATAACATCCATTTGTTTTTCAAGTTTTCCCATTAAGGTTCTTCAAAGCTACTTAAACATTTGTGGTTCTAGCATCC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:90 nM:i:4
chr6:59352460-59426290(-)_2578_3090_1:0:0_2:0:0_2d3     0       6       59352471        255     100M    *       0       0       ATTTTATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAAGATTTGTGTTTCTAGCATCCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:94 nM:i:2
chr6:59352460-59426290(-)_2624_3090_1:0:0_2:0:0_8b5     0       6       59352471        255     100M    *       0       0       ATTTTATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTGACATTAAGGTTCTTAAAAGGTACTTAAACATTTGTGGTTCTAGCATCCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:94 nM:i:2
chr6:59352460-59426290(-)_2644_3089_3:0:0_3:0:0_52b     0       6       59352472        255     100M    *       0       0       TTTTATCTTATTTGTCCCTGGGCTAACATACATTAGATTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:92 nM:i:3
chr6:59352460-59426290(-)_2491_3089_0:0:0_3:0:0_45c     0       6       59352472        255     99M1S   *       0       0       TTTTATCTTATTTGTCCCTGGGATCACATACATTTGTTTTTCAAGTTTTCACATTAAGGTTGTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTG    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:93 nM:i:2
chr6:59352460-59426290(-)_2578_3089_1:0:0_3:0:0_432     0       6       59352472        255     100M    *       0       0       TTTAATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACAATAAGGTTCTTACAAGCTACTTAAACATTTGTGGTTCTAGCATCCTC    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:92 nM:i:3
chr6:59352460-59426290(-)_2577_3088_0:0:0_1:0:0_648     0       6       59352473        255     100M    *       0       0       TTTATCTTATTTGTCCCTGGGATAACATACATTTTTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTCT    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:96 nM:i:1
chr6:59352460-59426290(-)_2628_3086_4:0:0_2:0:0_982     0       6       59352475        255     100M    *       0       0       TATGTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATAAAGGTTCTTAAAAGCTACTTAAACATTTGTGGTTCTAGCATCCTCTTG    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:94 nM:i:2
chr6:59352460-59426290(-)_2716_3086_1:0:0_3:0:0_8bc     1024    6       59352475        255     100M    *       0       0       TATCTTATTTGTCCCTGGGATAACATACATTTGTTTTTCAAGTTTTCACATTAAGGTTCTTAAAAGCTACTTAAACATTTGAGGTTGTAGCATCCACTTG    2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222    NH:i:1  HI:i:1  AS:i:92 nM:i:3

So only for one of the 3 groups of duplicates one of the reads (_8bc) is marked as a duplicate. Why are the other duplicates not regconised?

Your help would be much appreciated.

Many thanks in advance,

me

[Alexander Dobin]

Hi @ME,

testing your example I realized that I have a serious bug which prevents duplicate removal function from working with single-end alignments.

I will be fixing it, but it may take a while, as I need to make serious changes to the algorithm. In the meantime please use the samtools rmdup.

Cheers

Alex

[ME]

Dear Alex!

Thanks a lot for your reply!

I've tried a few other things in the meantime.

Samtools rmdup seems to miss duplicates if one of the two reads is soft-clipped. As far as I can see, it only considers the start coordinates of the alingments, not the true starts of the reads.

Picard MarkDuplicates is better in that respect, but has some other issues. It seems to compare the 5' portion of reads/alignments only, which can result in problems with gapped alignments. For instance, after using Picard on a test file similar to the one above (just much larger) there are still quite a few obvious duplicates not marked.

chr6:60731572-60829855(-)_599_1136_1:0:0_2:0:0_343      16      6       60732122        255     43M     *       0       0       GACATTCTTAGGCTTCAGGCTCATAGTCTTGGTAGCCTTCCTA     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:0  AS:i:42
chr6:60731572-60829855(-)_599_1061_1:0:0_1:0:0_78       16      6       60732122        255     40M942N3M       *       0       0       GACATTCTTAGCCTTCAGGCTCATAGTCTTGGTAGCCTTCCTC     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:1  AS:i:41


chr6:60731572-60829855(-)_596_1053_0:0:0_3:0:0_421      16      6       60732125        255     37M942N6M       *       0       0       ATTCTTAGGCTTCAGGCTCATAGTCTTGGTAGCCTTCCTCTGA     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:0  AS:i:43
chr6:60731572-60829004(-)_692_1197_2:0:0_0:0:0_e3       16      6       60732125        255     40M3S   *       0       0       ATTCTTAGGCTTCTGGCTCATAGTCTTGGTAGCCTTCCTATGA     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:1  AS:i:37


chr6:60731572-60829004(-)_691_1217_0:0:0_2:0:0_3da      16      6       60732126        255     36M942N7M       *       0       0       TTCTTAGGCTTCAGGCTCATAGTCTTGGTAGCCTTCCTCTGAA     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:0  AS:i:43
chr6:60731572-60829004(-)_691_1145_1:0:0_1:0:0_5c9      16      6       60732126        255     38M5S   *       0       0       TTCTTAGGCTTCAGGCTCATAGTCTTGCTAGCCTTCCTCTGAA     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:1  AS:i:35


chr6:60731572-60829855(-)_594_1116_2:0:0_2:0:0_35c      16      6       60732127        255     37M6S   *       0       0       TCTTAGGCTGCAGGCTCATAGTCTTGGTTGCCTTCCTCTGAAG     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:2  AS:i:32
chr6:60731572-60829004(-)_690_1166_1:0:0_2:0:0_2af      16      6       60732127        255     35M942N8M       *       0       0       TCTTAGGCTTCAGGCTCATAGTCTTCGTAGCCTTCCTCTGAAG     2222222222222222222222222222222222222222222     PG:Z:MarkDuplicates     NH:i:1  HI:i:1  nM:i:1  AS:i:41

Note that all these alignments are on the reverse strand and that in each pair one alignment is gapped. It seems, Picard only compares the 5' ends (i.e. 3' ends in genome orientation as reads are on the reverse strand), which are mapped to different regions, and hence fails to recognise these pairs as duplicates. Yet, the (major) 3' parts of these alignments are identical, and so they could be identified as duplicates quite easily. As a work-around, I've tried swapping forward and reverse flags after a first round of MarkDuplicates and then do a second round. This seems to work, but is quite messy of course.

So, it would be great to have a tool that can handle all these different situations for single-end alignments properly.

Thanks a lot for your help,

ME