Alternatives to Picard Mark Duplicates

[dpryan]

If he hasn't written one yet or posted a link to one written by someone else, then there's probably not been one written for gsnap's output (looking at the format, it'd be a pain to write one).

[fongchun]

Right. That's why I was wondering if there were alternatives to Picard mark duplicates that can handle these chimeric read formats.

I know that Thomas Wu, Gsnap author, said he was working on a package that could handle it, but I don't think it has been made public unfortunately...

[dpryan]

Well, picard can mark the duplicates, it just can't handle the chimeric read format of gsnap. So, you have some decisions to make...

[fongchun]

That's interesting. I've also read that it "depends" on your biological question especially if you are looking for SNVs. I would think that it is important to "mark" them and not necessarily remove them so that you can specify to the tool ahead of time how to handle it.

You are right in that a lot of tools break (including Picard). So it would seem that there is no way to mark duplicates in RNA-seq?

[dpryan]

It's typically not necessary to remove duplicates in RNAseq. With highly expressed genes, one would expect these.

BTW, a lot of tools will break when fed chimeric reads. More so since the SAM spec was only recently updated to include them and most mappers don't yet follow the spec.

[fongchun]

Hi,

Does anyone know any alternatives to Picard mark duplicates that is capable of handling split-reads in RNA-seq data?

I have a gsnap bam file that has been aligned and a read that it split is represented on two different lines (as indicates by a 'XT' flag. This ends up throwing an error when using Picard mark duplicates:

Exception in thread "main" net.sf.picard.PicardException: Value was put into PairInfoMap more than once. 6: test:M00569:47:000000000-A50YA:1:1107:13835:21068

The culprit lines in the bam file:

M00569:47:000000000-A50YA:1:1107:13835:21068    65      chr2    133036616       39      32H119M chr7    99156283        0       CTGCGGTTCCTCTCGTACTGAGCAGGATTACCATGGCAACAACACATCATCAGTAGGGTAAAACTAACCTGTCTCACGACGGTCTAAACCCAGCTCACGTTCCCTATTAGTGGGTGAAC      GHGFEEEEFFFHEHG2EGHGGCAD3FAGFDFGGFFEG3GFHHFFFHHGGHFHHFHFFF?3FGGGHFHHFHFDGFFHHFGGGGGFGGFFHFGGEHGHHGAGHHBHGHHHFGHF0GEDHHF      RG:Z:HEKRCOR1_clone1    MD:Z:T3CA8A1G17T11CATCAC1G23T3A6T9A1G18 NH:i:1  HI:i:1  NM:i:18 XW:i:6  XV:i:12 SM:i:39 XQ:i:40 X2:i:0  XS:A:-  XT:Z:GT-AG,0.87,0.98

M00569:47:000000000-A50YA:1:1107:13835:21068    81      chr4    76807259        39      119H32M chr7    99156283        0       GTTCAGACATTTGGTGTATGTGCTTGGCTGAG        GFFFD6GFDAGGGGGGFFGG4AAF?FFBBBBB     RG:Z:HEKRCOR1_clone1    MD:Z:11C17A2    NH:i:1  HI:i:1  NM:i:2  XW:i:0  XV:i:2  SM:i:39 XQ:i:40 X2:i:0  XS:A:+  XT:Z:GT-AG,0.87,0.98

M00569:47:000000000-A50YA:1:1107:13835:21068    129     chr7    99156283        38      119M32S chr2    133036616       0       CGTCGCCGGCAGTGCCACCGAGAAGCGCCGGCCTCGGGGCTGCCTACAGCGGCCCGGGAGAGGCTGTGGTGGGCCCGCGCGCGCGTGCGTAGGTGACAGGACACCGGCCGGGCCCGCCCTGGATAACTGGCTTGGGGCGGCCAAGCGTTCC      A?AAADAD10>>GGBBBF0EC0A/F1/AA/AEE/0E//>>/>/>1B0BBF@EG/E@/</</?/B/00BBFFC//?<////<---<-<.;E.-.<0CC00...9..;9-A-;@->--;@---AFF/B////;B---;@-@----99----9-      RG:Z:HEKRCOR1_clone1    MD:Z:9A32T29C4A5A3T15G5A9       NH:i:1  HI:i:1  NM:i:8  XW:i:8  XV:i:0  SM:i:38      XQ:i:40 X2:i:0  PG:Z:T

How do you typically handle duplicates in RNA-seq? Do you just use samtools rmdup?

Thanks,