3
$\begingroup$

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.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.