samtools piping problem

[malbuquerque]

Hello SEQanswers,

My ultimate goal is to pipe multiple different regions of a bam file to different commands in a single command line argument.

Basically, I am trying to generate some metrics from a BAM file that will be flagged for optical duplicates and then (without storing to disc) I want to pipe each chromosome to a program.

ex.

java -jar picard/MarkDuplicates.jar I=test.bam O=/dev/stdout | tee >(samtools view 'chr1' -| ./generate_metrics) >(samtools view 'chr2' -| ./generate_metrics) ...

any advice if my current implementation isn't the best?
The issue right now is attempting to pull a segment from the bamfile after an initial pipe.

Thanks,

Marco

[dpryan]

Well, that works, put the command is going to get really unwieldy. Why don't you either pipe to a perl/python/whatever script and have it pipe accordingly (the input to MarkDuplicates is already sorted, so your script would just need to pipe to one instance of ./generate_metrics at a time) or just rewrite ./generate_metrics to handle multiple chromosomes (assuming you have the source code).

[malbuquerque]

The idea is to parrallel these operations. So I will ultimately be submitting qsubs for each chromosomes operations, Its just a matter of segmenting the BAM file into its appropriate segments.

[dpryan]

I'd be a little surprised if piping like that into qsub would actually work. If you want to parcel this out to multiple nodes, why don't you just markDuplicates on a chromosome per node and generate_metrics on that?

[malbuquerque]

That sound like a great idea

Thanks! The duplicate algorithm (im assuming) is only looking for reads of the same length and position correct (i.e. there wont be duplicates accross different chromosomes).

[dpryan]

I think it just looks at the 5' position of both mates in a pair and marks duplicates if there are multiple pairs with the same 5' position(s). I could be misremembering that, though.

[malbuquerque]

So this is what I'm thinking from what you have told me

samtools view sample.bam 'chr1' | java -jar $PICARDROOT/MarkDuplicates.jar I= - O=sample_chr1.bam M=metrics_chr1.txt

Or is there something else I need to pass to I...

[dpryan]

Yeah, but I'd forgotten that MarkDuplicates can't take input via a pipe, so that won't work. I'm not coming up with a way to do that that doesn't end up involving considerable programming with MPI. Do you have the code for generate_metrics? It would seem easier to just write it to handle multiple contigs. Alternatively, presuming that actually marking the duplicates is the rate limiting step, then just run it on a single node.

[dpryan]

Oh, have a look at Biobambam, which can MarkDuplicates but writes its temp files to a local disk (it's geared more toward clusters). BTW, see also this thread on Biostars.

[rhinoceros]

Maybe gnuparallel could solve your parallelization? This thread explains it comprehensively..