I have a bam file with aligned reads from multiple experiments. For each experiment I have a portion of the sequence identifier, thus for example

  • HWI-D00310:448:CCG0KANXX:7:1101:2570:1976,

and is different from other the other of single-run flow cells in the analysis.

I want to demultiplex the input bam file into different bam files depending on the experiment. My only concern is a fast methodology to demultiplex a bam using combinations of single-run, flow-cell, and lane. What I came up with (the fastest) is something like:

# Parallelized per chromosome. Assuming only 2 plates here.
printf '%s\n' $(seq 1 22) X Y MT \
| parallel 'samtools view input.bam {} \
  | tee >(grep -f patterns_experiment1.txt > experiment1_chr{}.sam.noheader) \
        >(grep -f patterns_experiment2.txt > experiment2_chr{}.sam.noheader)' \
&& seq 1 2 \
| parallel 'samtools view -H input.bam && cat experiment{}_*sam.noheader \
  | samtools view -1 - > input_experiment{}.bam'

Is there a better/faster way to do it?

If I understood correctly read groups should be specified prior to mapping? I started from the already mapped BAM unfortunately, but looking at the original FASTQs I noticed the pattern in which some parts of the sequence identifier were unique for one experiment vs the others.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.