Samtools rmdup in a stream

[jdrum00]

I'm trying the following strategy to remove duplicates from small subsets of interest:

samtools view -b test.bam chr1:7660315-7660315 | samtools rmdup - - | samtools view -

...but it doesn't seem to remove anything. Am I doing something wrong syntactically, or do I perhaps have some of that data that rmdup doesn't work on? How do I tell if I have FR orientation, or if ISIZE is correctly set? (I'm new to using BAMs, and I don't create them myself.)

Here's a sample line.

938_1080_925 16 chr1 590592 0 30M * 0 0 TCATGTCAACTGCAAACAGAAACAATTTAA :BIIIICIIIIIIIIIHIIIIIIIIIIIII RG:Z:0 CS:Z:T003003011002211001312101211312 CQ:Z:?0>59;7=/;=><54<5<39?/@/5>3A):

(And yes, according to other posts, "I should really use Picard for this". But that has weird Java problems on my system that I don't have the time or permissions to fix.)

Thanks!

[lh3]

For samtools, you must set ISIZE correctly; otherwise it won't work. If you have single-end reads, use rmdup -s. To use Picard, you must have MRNM and MPOS set correctly for paired end reads (so far as I know). Picard is more theoretically correct about duplicates. You may see tiny differences between samtools and Picard results, but in practice this does not matter too much.

[jdrum00]

Thanks, but I don't know if I have ISIZE set correctly or if I have single-end reads. Those are my actual questions.

[Pepe]

You need to read the sam format specifications:
Look it up on google.
Column number 9 is the Inferred insert SIZE (ISIZE). If all of those in your file are 0, then you have single ends (or no read mapped to the reference )

[jdrum00]

You need to read the sam format specifications:
Look it up on google.
Column number 9 is the Inferred insert SIZE (ISIZE). If all of those in your file are 0, then you have single ends (or no read mapped to the reference )

I've read the spec backwards and forwards; it's a matter of understanding the details. I don't see anywhere that it says "if you have all zeroes, that means single ends, so use rmdup -s" -- so, thanks.

[lh3]

To tell if reads are paired in sequencing or mapped in pairs, one should look at the FLAG (2nd) column.

[jdrum00]

I know about the flags field, but I don't know what those things mean in detail. I've Googled for resources about pairing, but I don't think I'm getting the whole picture. I just don't know enough about sequencing technology yet, I guess.

Is the combination of pipes I'm trying reasonable, if I get the right options set up? That is, am I correct that I want a -b on my first view command, so that I'm sending BAM format to rmdup? And should those stdin/stdout dashes work, the way I have them? Thanks.